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ABSTRACT 


(Distribution  Limitation  Statement  No.  3) 

The  Modified  Elemental  Volume  Dose  Program  (MEVDP)  generates 
ordered  path-length  areal  densities  for  primary  electron,  electron- 
bremsstrahlung,  and  secondary  particle  radiation  transport  calculations. 

The  code  also  generates  standard-material  areal-density  distribution  functions 
for  proton  and  heavy  ionizing  nuclear  radiation.  The  primary  and  secondary 
areal-density  functions  can  be  used  for  particle  transport  calculations  to 
compute  emergent  fluxes  and  energy  deposition.  The  code  has  been  success¬ 
fully  run  with  the  complex  Apollo  command  and  service  modules  and  the  lunar 
module,  which  are  represented  by  1000  elemental  volume  shield  configura¬ 
tions.  The  MEVDP  has  evolved  as  a  versatile,  accurate,  and  fast  executing 
program. 

A  user's  manual,  which  is  an  integral  part  of  this  report,  illustrates 
the  procedure  for  input  data  preparation  and  execution  of  the  MEVDP. 
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SECTION  I 


INTRODUCTION 


To  facilitate  the  calculation  of  primary  radiation  transport  through 
the  complex  geometry  of  the  Apollo  vehicle,  the  North  American  Rockwell 
Corporation  Space  Division  (NR/SD)  developed  the  Elemental  Volume  Dose 
Program  (EVDP).  The  EVDP  was  developed  to  calculate  primary  radiation 
doses  from  protons  and  alpha  particles  associated  with  solar-flare  particle 
events  and  the  earth's  trapped  radiation  for  deep  space  and  earth  orbital 
missions.  The  EVDP  converts  the  material  paths  traversed  to  an  equivalent 
standard-material  path  length.  Because  of  the  general  versatility  and  com¬ 
putational  efficiency  of  the  EVDP,  a  modified  version  of  the  EVDP  has  been 
developed  for  electron,  electron-bremsstrahlung,  and  secondary  particle 
transport  calculations. 

The  modified  EVDP  (MEVDP)  generates  spatially  oriented  arrays  of 
material  type  and  associated  thicknesses  which  are  encountered  by  radiation 
traveling  toward  a  dosimeter  point  within  any  complex  geometrical  shielding 
configuration.  These  thickness  arrays  are  available  in  the  order  in  which 
the  incident  rays  encounter  the  different  materials.  The  modified  program 
also  retains  the  EVDP  capability  for  computing  the  standard-material  areal- 
density  distribution  function  versus  fractional  solid  angle  for  heavy  charged 
particle  dose  computations.  The  new  version  of  the  EVDP  contains  the 
source  ray  selection  and  geometrical  options  of  the  original  EVDP,  including 
the  composite  shield  routine.  These  features  have  been  modified  and  aug¬ 
mented  by  a  subroutine  to  order  the  computed  nuclear  transport  parameters. 
The  ordered  arrays  can  be  used  for  electron  and  secondary  transport  calcu¬ 
lations  with  the  straight-ahead  approximation. 

This  report  documents  the  following  aspects  of  the  modified  EVDP: 

(1)  geometrical  sectoring  technique,  (2)  analytical  techniques  incorporated 
in  the  code,  (3)  user's  manual,  (4)  modified  program  listing,  and  (5)  an 
AFWL  sample  problem  and  its  computer  solution.  The  analytical  techniques 
and  logical  operations  used  in  the  code  are  discussed  in  detail.  Analytical 
equations  are  derived  and  correlated  with  the  source  program  FORTRAN 
listing.  The  user's  manual,  an  integral  part  of  the  report,  concisely  and 
systematically  illustrates  the  procedure  for  data  preparation  and  its  execu¬ 
tion  to  obtain  computer  solutions.  The  solution  for  a  sample  problem 
illustrates  the  format  of  the  MEVDP  input  data  and  demonstrates  that  the 
program  is  operational. 
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SECTION  II 


GEOMETRICAL  SECTORING  TECHNIQUE 


The  dosimeter  position  can  be  imbedded  inside  any  of  the  system 
structural  components  or  the  spacecraft  crew.  The  source  flux  is  assumed 
to  be  isotropic  and  that  each  selected  particle  ray  maintains  its  original 
direction  when  penetrating  the  material  between  the  incident  particle  and 
the  dosimeter.  Both  random  and  systematic  methods  can  be  used  to  select 
the  ray  directions  of  the  source  nucleons  or  photons. 

Versatile  elemental  volumes  will  be  utilized  to  permit  a  relatively 
complete  description  of  the  spacecraft  and  crew  structural  configurations 
and  compositions.  This  versatility  will  eliminate  the  necessity  for  exces¬ 
sive  use  of  homogenization  in  the  geometrical  representation  of  subsystem 
components.  This  section  presents  the  technical  analysis  and  Section  III 
directly  correlates  the  analysis  to  the  computer  program.  The  analysis 
begins  by  presenting  the  method  for  selecting  the  directions  of  the  source 
nucleons  or  photons. 

1.  SELECTION  OF  THE  INCIDENT  SOURCE  DIRECTION  COSINES 

If  the  random  method  is  used  for  dose  or  flux  computations,  the 
direction  cosines  will  be  selected  by  random  sampling  uniformly  over  the 
unit  sphere.  The  systematic  method  will  select  directions  at  the  center  of 
a  preselected  number  of  equivalent  solid  angles.  Also  different  solid-angle 
increments  can  be  used  for  ray  tracing  in  different  portions  of  the  unit 
sphere  centered  at  the  dosimeter  position.  The  differential  selection  of 
the  solid-angle  increment  can  result  in  a  more  accurate  dose  computation 
by  selecting  larger  ray  densities  for  the  higher  dose  regions. 

a„  Random  Direction  Cosines 


For  random  selection  of  unnormalized  ray  direction  cosines  (oQ,  pD, 
Yq)  over  the  unit  sphere 


where 


\  =  1  ■  2t>1 


\  -  1  - 


-  1  -  2r), 


</  +  Pn  +  V2  <  1 

o  “o  1  n 


2 


and  Hi»  rj 2 »  and  H  3  are  random  numbers  from  a  uniform  distribution  in  the 
interval 


0  <  q.  <  1  ;  i=  1 ,  3 
1 

The  normalized  direction  cosines  are  a,  p,  and  y  where 

a 

P 

Y 

€ 

b.  Systematic  Selection  of  Direction  Cosines 


(1) 


The  ray  direction  cosines  are  selected  at  the  center  of  equal  solid- 
angle  increments  as  depicted  in  Figure  1.  A  right-handed  spherical  coor¬ 
dinate  system  is  used  with  6  and  <f>  as  colatitude  and  azimuthal  angles, 
respectively. 

A  microscopic  solid-angle  portion  of  the  unit  sphere  is  considered, 
which  is  defined  by  the  colatitude  interval  (6 j,  0p)  and  the  azimuthal 
interval  (0j,  <p p)  and  divided  into  NSA  microscopic  solid-angle  increments. 
The  following  definitions  are  used: 

NTS  =  total  number  of  solid-angle  increments  in  the  unit  sphere 

FSA  =  fraction  of  the  total  unit  sphere  solid  angle  in  the  angular 
region  (Gj,  0F;  <pF) 

NSA  =  number  of  solid-angle  increments  in  the  FSA  portion  of  the 
unit  sphere 

Nq  =  number  of  0  increments  in  NSA 

N  =  number  of  <b  increments  in  NSA 

0 


3 


4 


For  symmetry,  the  6  and  <p  angular  magnitudes  of  all  incremental  solid 
angles  are  considered  to  be  equivalent,  or 

♦f  -*1  eF  -  6I 


The  fraction  solid  angle  (FSA)  can  be  derived  as 


FSA  = 


if'/ 


sin  0d0  d0 


♦i  ei 


(^r1)  (cos  ei  ■ 


cos  eF) 


NSA  =  FSA  •  NTS  =  N  N 

e  0 


From  Equations  2  and  4 


n2  (**  '  01 

=  Ne  W  - 


N'  =  IPO 
0 


l(^r) 


NSA  +  1 


where  IPO(X)  is  defined  as  the  integral  part  of  (X).  The  modification  of 
Nq  (Equation  5) is  necessary,  because  Ng  must  be  a  nonzero  integer. 

Since  the  incremental  solid  angles  are  equivalent 


1*1  -  ^cos  e;  -  cos  eKj  ~  I^os  8; 


-  cos  eF\  /♦,  -  «U 


where  K  =  1,  2,  ....  Ng. 
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It  follows  that 


cos  0  =  cos  6  -  — r  (cos  0  -  cos  0_) 


'K 


I  N, 


6 


I 


1 


eK  =  cos  (coseK) 


(7) 


The  6j^  are  the  6  boundaries  of  the  incremental  solid  angles. 
From  Equation  2 


N'  =  IPO 
<P 


The  <p  boundaries  are  defined  by 


/*F  -  V 

\  F  -  «I, 


(Ne  -  i) 


+  l 


*|  =  V  <*F  '  *I>  ;  l=°’  Ni 


(8) 


The  solid-angle  midpoints  (0  s  0j^)  and  ( <b  2  are 


8MK  -  °-5|eK  +  (W 

♦M,  =  0.  5  +  <fi+1)  (9) 

The  total  number  of  computed  incremental  solid  angles  is  the  product  of  Nq 
and  N^. 

Finally,  the  systematic  direction  cosines  are 

a  =  sin  0  cos  <f> 

(}  =  sin  0  sin  <p  (10) 

Y  =  cos  0 

The  next  logical  step  is  utilization  of  the  ray  direction  cosines  to 
compute  the  traversal  path  lengths  through  the  geometrical  structure. 
Several  different  coordinate  systems  will  be  used  to  simplify  the  path-length 
computations. 
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2.  COORDINATE  SYSTEMS 


The  geometrical  structural  configurations  are  originally  defined  by 
elemental  volume  coordinates  and  parameters  in  a  coordinate  system 
(Figure  2)  designated  as  the  Absolute  Coordinate  System  (ABCS).  In  order 
to  substantially  decrease  the  amount  of  ray  tracing  or  tracking  required, 
the  elemental  volumes  are  classified  according  to  their  octant  location  in 
a  coordinate  system, centered  at  the  detector. 

With  this  choice  of  the  coordinate  origin,  the  incident  particle  track 
must  traverse  only  one  octant. 

However,  if  the  volume  octants  are  selected  relative  to  the  ABCS, 
there  will  be  unnecessary  tracking  because  the  source  rays  are  directed 
toward  the  dosimeter,  which  is  generally  displaced  from  the  origin  of  the 
ABCS.  Therefore,  the  proton  track  could  traverse  four  octants.  Tracks 
along  the  coordinate  axes  are  considered  as  special  cases. 


Figure  2.  Absolute  Coordinate  System  (ABCS) 
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The  pertinent  parameters  and  coordinates  used  to  define  the  elemental 
shield  volumes  are  depicted  in  Figure  3  for  the  ABCS. 

Figure  4  depicts  the  Dosimeter  Coordinate  System  (DSCS)  and  the 
ABCS.  To  transform  a  coordinate  point  Ra  =(xA»  YA*  za)  °f  *he  spacecraft 
to  the  DSCS,  the  following  displacement  is  required: 

rd  =  ra  -  rad 

XD  =  XA  -  XAD 

(ID 

YD  =  VA  -  YAD 

ZD  =  ZA  '  ZAD 

where 

R^  =  (xa,  YA*  ZA)  =  position  of  the  spacecraft  point  relative 
to  the  ABCS 

Rd  =  (XD»  Yd»  zd)  =  posifio11  °f  the  spacecraft  point  relative 
to  the  DSCS 

Rad  =  (XAD»  YAD*  ZAD)  =  position  of  the  DSCS  origin  relative 
to  the  ABCS 

For  several  elemental  volumes,  it  is  most  convenient  to  have  a 
detector  coordinate  system  with  one  or  more  of  its  axes  directed  along  the 
reference  axes  of  the  elemental  volumes.  Since  the  new  detector  system 
has  the  same  origin  as  the  DSCS,  the  direction  cosines  can  be  converted  to 
the  new  system  by  matrix  rotations.  The  new  system  will  be  designated  as 
the  Rotated  Detector  Coordinate  System  (RDCS).  The  RDCS  is  illustrated 
in  Figure  5  for  a  cylinder.  6  and  4>  are  the  component  rotation  angles.  The 
4>  rotation  about  the  DSCS  z  axis  is  followed  bya0  rotation  about  intermediate 
coordinate  system  (xj,  yj,  Zj)  negative  Xj  axis. 

In  order  to  select  the  rotation  matrix  transformations  for  the  elemental 
volumes,  it  is  necessary  to  determine  convenient  coordinate  systems  relative 
to  the  elemental  volume  reference  axes.  The  most  convenient  axes  were 
selected  as  those  which  facilitate  calculations  of  path  lengths  to  the  elemen¬ 
tal  surfaces. 
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Dosimeter  (DSCS)  and  Absolute  (ABCS)  Coordinate  Systems 


3.  SURFACE  PATH-LENGTH  COMPUTATIONS  FOR  PLANES  AND 
CONIC  SURFACES 

Path  length  computations  are  described  for  the  following  elemental 
surfaces:  sphere,  ellipsoid,  cone,  cylinder,  and  plane. 

The  path  length  is  generally  computed  as  the  difference  between  the 
two  path  lengths  to  the  surfaces  of  the  volumes.  If  the  dosimeter  is  inside 
the  volume,  the  volume  path  length  is  determined  by  a  single  path  length. 
However,  additional  tests  are  required  for  truncated  volumes  and  for  the 
hexahedron:  these  are  discussed  in  the  latter  portion  of  this  section. 

The  equation  of  the  surfaces  which  encompass  the  elemental  volumes 
are  quadratic  forms.  The  hexahedron  and  truncated  surfaces  are  combina¬ 
tions  of  planes.  Figure  6  shows  the  points  and  reference  axes  used  to  define 
the  ellipsoid.  This  volume  is  used  to  illustrate  the  method  of  computing 
path  lengths  for  all  volumes  except  the  hexahedron. 

The  following  definitions  are  used  in  Equations  12  through  17  (see 
Figure  6): 

(xr,  Yr,  Zr)  s  coordinates  (RDCS)  of  the  point  where  the  particle 
path  traverses  the  ellipsoid  surface 

(XR4,  yR4,  zr4)  2  center  of  the  ellipsoid 
(a,  b,  c)  s  ellipsoid  major  axes 
(x0.  yQ,  zQ)  z  a  point  on  the  particle  path 
(or,  Pr,  Yr)  *  direction  cosines  of  the  particle  path  in  the  RDCS 

In  the  rotated  coordinate  system,  the  equation  of  the  ellipsoid  surface 
is 


The  equation  of  a  line  through  the  point  (xD,  yQ,  zD)  in  the  RDCS  is 


12 


where 


Since  the  path  length  to  the  ellipsoid  from  the  point  (xQ,  y0,  zQ)  is  P, 


P  = 


/z  ^  -  z 
'  R  o 


rR 


(14) 


Equation  13  can  be  used  to  derive 


|zR  '  zo> 
XR  =  ’to+“R-^— 


(ZR  -  V 


yR  =  yo  +  ?R  V[ 


ZR  Zo+VR  V 


,ZR  '  V 


When  Equation  14  is  substituted  into  Equation  15 


(15) 


*R  =  xo  +  “RP 
yR  =  yo  +  PRP 
ZR  =  zo+yRP 


(16) 


and,  subsequently,  when  Equation  16  is  substituted  into  Equation  12,  it  follows 
that 

2 

AP  +  BP  +  C  =  0 


where 


2  „  2 

fP 


-G) '(?)  ♦£) 


B  =  2 


°R  (Xo  ~  V  .  *  R  <yo  '  yR4>  .  VR  ,Zo  '  ZR4> 

2  +  2  +  2 
a  b  c 


(17) 
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and  the  two  path  lengths  are 


-B  ±  yB  -  4 AC 
2A 


The  solutions  for  the  cone,  cylinder,  and  sphere  are  also  quadratic 
in  form.  In  some  cases,  it  is  possible  that  A  =  0  or  one  path  length  is  an 
indeterminate  form.  For  example 


+  0 


Using  L'Hopital's  rule,  this  indeterminate  form  can  be  evaluated  as 


rojki±^!-4Ac 

Au-mo  h  H 


For  the  other  volumes,  A,  B,  and  C  are  as  follows: 


a.  Sphere 


A  =  1.  0 


B  ’  -2  taD*Dl+,5DSrDl+VDZDl) 


C  =  x_.  +  y_.  +  z_,  -  R 
Dl  7D1  Dl 


R  =  (XD2  -  XD1 


)  +(yD2-yDl)  +(ZD2'ZDlj 


b.  Cylinder 


a  =<*l  +az 

R  R 


B  =  2[a'R(Xo-XRl)+BR(''o-VRll] 


C  =  /x  -  x 


o‘XRl)  +(yo’yRI 


/  \2  T,2 

-  y^i  -  R 
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c.  Cone 


2  „2  2  2ft 
A  <>R+|iR-VR,in  6 


B  -  2  U  (x  -  *D  J  +  (y  •  yDJ  •  Vr,  <z  -  z  )  x  tan  0 
1  R  o  R2  R  o  R2  R  o  R2 


C  = 


(*o  -  XR2)‘  +  (yo  '  VRzf  -  iZo  -  ZR2 f  ta"2e 


(22) 


where  subscripts  R  and  D  refer  to  the  RDCS  and  DSCS,  respectively,  and  the 
numbers  used  in  the  subscripts  are  defined  in  Figures  7  and  8.  If  the  ray 
tracing  direction  passes  through  the  origin  of  the  coordinate  system 

x  =  y  =  z  =0 
o  o  o 


d.  Plane 


The  plane  path-length  equation  is  not  a  quadratic  function.  In  Fig¬ 
ure  9,  a  plane  defined  by  a  point  P  and  a  unit  vector  "epL  is  intersected  by 
a  line  which  is  defined  by  the  point  0  and  the  unit  vector  e^.  The  path  length. 


which  is 


TH 


s  PTu  can  be  computed  from  the  equivalents  of  two  scalars 


( Rp 


-Ro'' 


p  r  P  p  *  p 

PL  TH  L  PL 


TH 


PL 


(VR 


o) 


e  •  p 

L  PL 


(23) 


where  PJH  =  0  ife^  epL  S  0 


(24) 


e.  Hexahedron 

The  hexahedron  is  a  combination  of  planes,  so  the  plane  path-length 
formula  can  be  used  to  determine  the  hexahedron  path  lengths^  The  hexa¬ 
hedron  plane  in  Figure  10  is  defined  by  vectors  V(l)  through  V(4).  Fig¬ 
ure  11  shows  a  projection  of  a  hexahedron  surface  in  the  (x^,  yp)  plane  of 
the  RDCS.  The  RDCS  is  defined  with  the  tracking  ray  directed  along  the  zp 
axis.  Points  B  and  A  represent  hits  or  misses  of  the  projected  hexahedron 
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Figure  9  .  Plane  and  Line  Vectors 
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face  for  the  origin  of  the  RDCS.  If  the  origin  is  inside  the  projected-plane 
perimeter,  the  ray  hits  the  face.  Figure  II  illustrates  a  set  of  angles  (6) 
which  will  be  used  as  a  test  to  determine  if  the  ray  hits  the  hexahedron  face. 
A  hit-or-miss  criterion  can  be  derived  from  the  sum  of  the  angles  sub¬ 
tended  by  the  sides  of  the  projected  plane  as  follows: 


4 

^  6^  =  0  ;  Miss 

i  =  l 

=  2tt  radians  ;  Hit  (23) 

However,  a  simpler 
simpler  hit-or-miss 

4 

I 

i=  1 


where 

V5>  ■  V11 

sign  (x)  «  |^| 


criterion  is  possible  because 
criteria  is 

<  180  deg] 

sign  * 

V‘» » Vi:'"|) 

*  4  ; 

Miss 

=  4  ; 

Hit 

(26) 


This  criterion  is  used  to  determine  which  faces  of  the  hexahedron  are 
intersected  and  Equation  23  is  used  to  compute  the  associated  path  lengths. 

4.  COORDINATE  SYSTEMS  FOR  PATH-LENGTH  COMPUTATIONS 

The  convenient  coordinate  systems,  and  directions  of  the  RDCS  z-axis 
orientations  are  defined  in  Table  I. 

Several  of  the  elemental  volume  types  require  transformations  from 
the  Absolute  Coordinate  System  (ABCS)  to  the  Rotated  Detector  Coordinate 
System  (RDCS). 
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rPT*  1  * V  T«  ^ 


Table  I 

Coordinate  Systems 


Elemental  Volume 

Coordinate 

System 

Direction  of  the  Axes 
Relative  to  Elemental  Volume 

Sphere 

DSCS 

- 

Cylinder 

RDCS 

Axis  of  cylinder  along  z 

R 

Cone  and  trun- 
c'  tod  cone 

RDCS 

Axis  of  cone  along  z^ 

Elhpsoid  and 
"'•iiicated  ellipsoid 

RDCS 

Axis  of  z^  along  truncation 
and  all  RDCS  axes  along 
major  axes 

Hexahedron 

RDCS 

z^  along  direction  of 
tracking  ray 

Hemisphere 

RDCS 

z^  along  direction  of 
truncation 

5.  COORDINATE  SYSTEMS  ROTATION  TRANSFORMATIONS 

Several  elemental  volumes  have  surfaces  for  which  path  lengths  can 
be  more  conveniently  computed  in  the  RDCS  by  transformations  which 
require  equivalent  matrix  rotations  of  the  DSCS  to  the  RDCS.  This  sub¬ 
section  derives  the  required  transformations. 

a.  Ellipsoid 


Since  the  path  lengths  are  computed  in  the  Rotated  Detector  Coordinate 
System  (RDCS),  the  ellipsoid  detector  system  coordinates  must  be  trans¬ 
formed  to  the  RDCS  as  follows: 

Let  R^  be  a  vector  equivalent  of  the  kth  point  of  the  ellipsoid  (Figure  6) 

and 


Rk  2R(k)  ss[x(k)lf  x(k)2,  x(k)3J  (27) 

Using  the  summation  convention  definition 

3 

AiBi 

i=l 
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it  follows  that  since  an  actual  vector  or  tensor  of  rank  1  is  invariant  in  any 
coordinate  system,  for  the  DSCS  and  RDCS 


R{k)  =  x(k)Ro  =rDox(k)Da  (28) 

where  subscripts  R  and  D  refer  to  the  RDCS  and  DSCS,  respectively,  and  the 
subscript  o  refers  to  the  oth  coordinate  component. 

The  ith  component  of  the  kth  point  in  the  RDCS  is 


x(k)Ri  =  (*Ri  *  iDo)  x(k>] 


Since  the  ellipsoid  points  numbered  k  =  1  to  4  in  Figure  6  can  be  transformed 
to  the  DSCS  by  Equation  11,  the  unit  vector  in  the  RDCS  can  be  resolved  in 
the  DSCS  as  follows  (Figure  12): 

Td  [x(k)D  -  x(4)d  ] 

lRk  (  3  ?  i  1/2  ;  k  =  e*  3 


3 

^  |x(k>Dj  -  x(4>Dj] 

j  =  l 


iRk  =  iDo  2(k)Da 


where  i  is  a  unit  vector  directed  from  point  4  to  point  k. 
Rk 


From  Equations  29  and  30 


r,k,Rj  *  19)3  'WD0  '  TE*1  X(k)Dff 


=  %  ,(i,D3  x(k)Eb  =  X(k,D« 


a=  X(1)D1  ’  X(4)D1  etC- 


b.  Hemisphere,  Cylinder,  and  Cone 

From  Figure  5  the  cone  coordinate  transformation  can  be  effected  by 
rotations  about  the  zp  and  xj  axes  of  angles  c  and  -6,  respectively.  If  a 


matrix  rotation  is  defined  by  )  where  refers  to  the  axis  system  and 

«  to  the  angle  of  rotation,  then  a  vectoi  position  5T  is  transformed  by 


x(RDCS)  =  MxJ(-e)  V  (0-90°)  x(DSCS) 


(32) 


and 


M  A(») 


cos  0  sin  0  0 

L -  sin  0  cos  0  0 

0  0  1 


M  (6) 

— xl 


The  x  matrix  is  a  column  matrix. 


0  0 
cos  6  sin  0 

-sin  0  cos  0 


Since 


cos  (0-90°)  =  sin  <p 
sin  (0-90°)  =  -  cos  <p , 


(33) 


M  ("6)  -zA  (*~9°0) 


sin4» 

cos  6  cos«J> 
sin  6  cos0 


-cos  4> 
cos0  sin0 
sin  6  sin  0 


e  M 


(34) 


The  trigonometric  functions  of  0  and  4>  are  determined  as  follows: 
If  the  direction  cosines  of  the  centerline  of  the  cone,  cylinder,  etc.  are 
defined  by  (a,  p,  y),  then 


cos  0  =  Y 

sin  6  =  ^1  -  Y^) 


since  0  £  180  degrees. 

If  o=0,  sin  0  =  1.  0  x  sign  of  (3  and  cos  0  =  0.  0. 
If  o=0,  and  Y  =  1»  sin  0  =  cos  0  and  cos  0  =  0.  0. 


(35) 
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Otherwise 


where 


tan  <p 
cos  <p 


(3  /  a 

2 

±(1  +  tan  <f>) 


/  2 


cos  0  >  0 
cos  <p  <  0 
sin  0  =  ± 


if  <*  >  0 

if  a  <  0 
2 

(1  -  cos  <p) 
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(36) 


where 


sin  0  <  0  if  (3  <  0 
sin  0  >  0  if  (3  >  0 

The  matrix  equivalent  of  Equation  34  is  used  to  transform  the  direction 
cosines  as  follows: 


=  MaD  (37) 


where  (M)jj  =  f(i)j}j  for  the  ellipsoid.  The  discussion  of  the  analysis 
required  to  compute  the  path  lengths  to  the  elemental  volume  surfaces  is 
now  complete.  However,  these  surface  path  lengths  must  be  transformed 
to  volume  path  lengths;  i.  e.  ,  path  lengths  through  volumes. 

6.  VOLUME  PATH-LENGTH  COMPUTATIONS 

The  elemental  volumes  are  in  general  composite  surfaces,  including 
conical  and  plane  surfaces.  It  is  necessary  to  know  whether  the  dosimeter 
is  inside  the  elemental  volume,  because  in  some  cases  the  path  through  the 
volume  is  the  difference  between  two  allowed  path  lengths,  while  if  the 
elemental  volume  contains  the  dosimeter,  the  path  through  the  volume  is 
equivalent  to  one  of  the  allowed  surface  path  lengths. 


Q 

-R 


a 


R 


Pi 


LYR 
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a.  Test  for  Dosimeter  Inside  Elemental  Volume 


Sphere 


The  dosimeter  is  inside  when  the  distance  from  the  dosimeter  to  the 
sphere  center  is  less  than  the  sphere  radius,  or  (Figure  7) 

XD(1)  +  yD(1)  +  ZD(1)<  ^  (38) 

This  condition  for  the  conic  surfaces  corresponds  to  C  less  than  zero 
(Equations  17  and  20  through  22). 

He  misphere 

A  sphere  with  the  same  radius  as  the  hemisphere  and  centered  at  the 
hemisphere  center  must  contain  the  dosimeter,  and  the  origin  components 
in  the  RDCS  must  be  bounded  by  the  z  components  of  points  1  and  2  (Fig¬ 
ure  8).  This  is  true  if 


zR(l)  zR  (2)  <  0 

and 

xr2(1)  +yR2(l)  +  zr2(1)  <  R2 


(39) 


Cone 


The  origin  must  be  bounded  by  the  z  components  of  the  base  and  apex 
of  the  cone,  and  the  circular  intersection  of  the  (xr,  yR)  plane  and  the  cone 
must  contain  the  dosimeter  (Figure  8)*  These  conditions  are  satisfied  if 

zR(l)  zR (2 )  <  0 


and 


xr(2)  +  yR(2)  <  zR(2)  tan  0 


(40) 


Similar  logic  can  also  be  used  to  determine  the  following  inequalities  for  the 
cylinder  and  ellipsoid. 

Cylinder  (Figure  8) 


ZRU)  zr(2)  <  0 


28 


and 


4(1>  +  Ypd)  <  R 


(41) 


Ellipsoid  (Figure  6) 


4(4)  4<4)  zr(4) 

___  +  -j-  +  — 

a  b  c 

and  for  double  truncation  also 


(42) 


*R(S)  zr(6)  <  0 


Hexahedron 


If  M  is  the  midpoint  of  a  vector  from  points  1  to  5,  then  the  dose  point 
D  is  inside  if  MD  is  less  than  MP  (Figure  13).  P  is  the  point  of  intersection 
of  a  vector  from  M  in  the  direction  of  D. 

After  it  is  determined  that  the  dosimeter  is  outside,  additional  logic  is 
necessary  to  accommodate  unique  cases  such  as  glancing  incidence,  etc. 

b.  Volume  Path-Length  Logic 


There  are  three  basic  types  of  computed  path  lengths,  i.  e. ,  positive, 
negative,  and  path  lengths  to  unacceptable  regions  of  a  surface.  Figure  14 
illustrates  these  path  lengths.  The  particle  ray  direction  is  defined  from 
the  dosimeter  position  inside  the  cone  to  point  A  on  the  cone  surface.  This 
ray  also  intersects  the  cone  base  plane.  The  distance  from  the  dosimeter 
to  point  A  (Pj)  is  defined  as  a  positive  path  length.  The  path  length  from  the 
dosimeter  to  the  base  plane  is  an  unacceptable  path  length  because  only  the 
base  plane  path  lengths  to  the  circular  region  of  the  base  plane  are  allowed. 
The  path  lengths  computed,  for  this  case,  with  Equation  18  will  have  positive 
and  negative  values.  The  negative  path-length  value  will  be  P2  in  Figure  14. 
This  is  another  case  of  an  unacceptable  path  length.  The  unacceptable  path 
lengths  are  used  in  the  path-length  logic  and  are  assigned  arbitrary  large 
values  of  10^®. 

If  the  dosimeter  is  outside  a  composite  surface  volume,  in  general, 
the  difference  between  two  path  lengths  is  the  path  length  through  the  volume. 
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Pj  ■  POSITIVE  PATH  LENGTH 
P2  «  NEGATIVE  PATH  LENGTH 


Figure  14.  Illustrated  Positive  and  Negative  Path  Lengths 


This  is  also  true  when  tho  tracking  ray  is  tangent  to  one  of  the  volume  sur¬ 
faces  and  the  path  lengths  are  equal.  However,  if  the  ray  passes  through  a 
single  surface  and  the  intersection  of  two  of  the  composite  surfaces,  there 
will  be  three  path  lengths  with  two  being  equal.  Therefore,  for  all  volumes 
except  the  hexahedron,  if  there  are  three  or  more  paths,  the  minimum  path 
is  determined;  then  it  and  all  equivalent  paths  are  set  equal  to  10^0,  This 
is  repeated  for  tho  second  minimum,  and  the  volume  path  length  is  the 
difference  between  the  two  minimum  paths. 

For  the  hexahedron,  if  there  are  only  two  allowed  path  lengths,  the 
volume  path-length  computation  is  the  same  as  for  the  other  elemental 
volumes.  If  there  are  three  or  more  allowed  path  lengths,  the  minimum  of 
the  path  lengths  is  picked,  and  it  and  all  equivalent  values  are  set  to  1020, 
Next,  the  second  minimum  is  picked  and  tested  for  equivalences  to  lO2^. 

If  the  second  minimum  is  equivalent  to  1020,  there  is  glancing  incidence,  and 
the  volume  path  length  is  zero.  Otherwise,  the  volume  path  is  the  difference 
between  the  first  and  second  minima. 

c.  Composite  Shield 


The  elemental  volumes  are  defined  as  positive  and  negative  or  solid 
and  void,  respectively.  A  negative  or  void  volume  is  defined  as  voiding  all 
the  portions  of  the  solid  volumes  which  intersect  its  spatial  region.  The 
negative  volume  only  voids  solid  shields  in  its  particular  composite  shield. 
The  composite  shield  is  a  combination  of  up  to  10  elemental  volumes. 

The  composite  shield  permits  imbedding  in  which  the  positive  shield 
may  contain  a  portion  or  all  of  the  negative  or  void  shield.  Also,  the  com¬ 
posite  shield  positive  shield  components  may  have  different  material  types 
and  densities.  Portions  of  several  of  the  positive  shield  components  may 
contain  portions  of  a  single  negative  or  void  shield.  These  additions  result 
in  a  very  significant  improvement  in  the  geometrical  representation  capa¬ 
bility  of  the  progiam. 

The  negative  or  void  shields  are  defined  such  that  the  portions  of 
positive  shields  which  are  contained  in  the  negative  shield  volumes  are 
considered  as  void  regions.  The  positive  shields  are  considered  separately 
and  their  portions  contained  in  the  negative  shields  voided.  Since  the  indivi¬ 
dual  negative  shields  may  intersect  some  or  all  of  the  positive  shields,  all 
the  negative  shields  must  be  tested  for  intersections  with  each  positive  shield. 
The  tracking  ray  will  sequentially  traverse  the  negative  shields  according  to 
their  increasing  distance  from  the  dosimeter  location.  So  the  negative  shields 
are  ordered  in  a  path-length  array  [PN  (I,  J)]  where 

I  refers  to  a  set  of  path  lengths  for  a  negative  shield 
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J  refers  to  the  two  path  lengths  at  the  tracking  ray’s  entrance  and 
exit  from  the  elemental  negative  shield,  i.  e.  ,  [PN  (I,  1)  PN  (1,2)] 

and 


PN  (I  +  1,  1)  >  PN  (I,  1),  1=1,  No.  negative  shields  traversed  in  the 
composite  shield 

The  positive  shield's  path  lengths  are  redefined  in  an  array  PP(I,  J)  where 
I,  J  refer  to  the  positive  shields  and  their  two  path  lengths,  respectively. 

The  array  of  positive  path  lengths  is  not  ordered.  Therefore,  they  are 
redefined  in  the  same  sequence  that  the  respective  positive  shields  are 
tracked.  Another  array  IMAT  (I)  contains  the  material  type  corresponding 
to  the  positive  shield.  IMAT  refers  to  the  material  type  for  the  areal- 
density  computions.  In  Figure  15,  the  logic  for  the  traversel  path  lengths 
and  plotting  points  is  illustrated  for  the  positive  shield  (+1)  and  the  void 
shields  (-2  and  -3).  For  convenience  the  positive  and  negative  shield  path 
lengths  are  defined  as  [P  (1),  P  (4)]  and[P  (2),  P  (3)],  respectively  for  the 
first  negative  shield  intersected.  If  the  negative  (-3)  shield  does  not  inter¬ 
sect  the  positive  shield  (+1),  the  following  test  must  be  satisfied: 

P  <3)  <  P  (1)  or  P  (2)  >  P  (4) 

Otherwise,  the  negative  voids  a  portion  of  the  positive  shield. 

If  P(2)  >  P  (1),  the  portion  of  the  positive  shield  between  P  (1)  and  P  (2)  is 
plotted  or  the  path-length  equivalent  areal  density  is  computed  and  added  to 
the  other  equivalent  path  lengths  of  the  particular  tracking  ray.  However, 
if  P  (2)  <  P  (1)  the  negative  shield  can  completely  void  the  positive.  In  this 
case,  if  P  (3)  >P  (4),  the  positive  shield  track  is  completely  voided  and 
another  positive  shield  is  selected.  However,  if  P  (4)  >  P  {3),  the  portion 
of  the  positive  shield  between  P  (1)  and  P  (3)  is  voided  by  redefining  the 
positive  shield  path  lengths  as  Q  (1),  Q  (4),  where 

Q  (1)  =  P  (3) 

Q  (4)  =  P  (4) 

If  P  (2)  >  P  (1),  the  test  must  also  be  made  with  P  (3)  and  P  (4)  to  determine 
if  the  remainder  of  positive  shield  is  voided.  If  all  the  negative  shields  are 
processed  without  satisfying  the  condition  P  (3)  2  P  (4),  Q  (3)  >  Q  (4)  etc. 
the  path  length  between  P  (3)  and  P  (4)  is  used  and  another  positive  shield 
is  selected. 
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Composite  Shield  Tracking  Sequence 


d.  Hit-or-Miss  Criteria 


To  eliminate  unnecessary  path -length  computations,  the  first  portion 
of  the  volume  path -length  calculations  will  determine  if  the  tracking  ray 
misses  a  simple  boundary  which  contains  the  elemental  volume.  For  example, 
the  sphere  contains  the  hemisphere  elemental  volume.  However,  if  the  ray 
hits  the  spherical  volume  it  may  hit  or  miss  the  hemisphere. 

Sphere 


The  sphere  is  missed  if  the  discriminant  of  Equation  18  is  negative 

B2  -  4AC  5  0 

which  means  the  path  length  through  the  sphere  is  imaginary  or  zero. 
Cylinder  and  Cone 


This  criterion  is  based  on  intersection  of  the  tracking  ray  with  the 
projection  of  the  cylinder  and  the  larger  radius  of  the  cone  on  the  (x,  y)  plane 
in  the  RDCS.  This  is  done  by  using  the  center  point  of  the  circular  projection 
as  the  center  of  a  sphere  with  the  radius  of  the  projection.  The  tracking  ray 
misses  the  volume  if  the  component  of  the  tracking  ray  in  the  (x,  y)  plane  does 
not  intersect  the  aforementioned  sphere.  The  discriminant  of  Equation  18 
becomes 


B1  •  +  »*,]/(•**♦  IV,1  ),y* 

_  2  2  2 
C1  *R1  +  yR  1  "  R 


and  there  is  a  miss  if 


B 


-  4Aic, < 


(43) 


Ellipsoid 


The  ellipsoid  miss  criterion  is  the  same  as  that  for  the  sphere,  except 
that  the  ellipsoid  is  in  the  RDCS  and  the  radius  is  the  largest  ellipsoid  major 
axis. 
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Hexahedron 


The  hexahedron  and  each  plane  face  of  the  hexahedron  is  projected  on 
a  plane  perpendicular  to  the  ray  direction.  Then,  in  the  RDCS,  a  rectangle 
is  defined  by  the  minimum  and  maximum  x  and  y  values  (Figure  16).  If  the 
rectangle  does  not  contain  the  origin,  the  tracking  ray  misses  the  face. 

Another  form  of  hit-or-miss  check  is  tracking  by  octant,  and  consider¬ 
ing  only  shields  in  the  octant  of  the  tracking  ray. 

The  octants  are  defined  by  the  DSCS  coordinate  sign  convention  in 
Table  II.  By  using  the  DSCS,  the  matrix  rotation  calculations  required  for 
the  RDCS  can  be  eliminated  if  the  ray  and  shield  octants  are  not  identical. 

Table  II 

Octant  Definition 


Coordinate  Component  Signs 

Octant  No. 

X 

y 

z 

- 

- 

+ 

1 

2 

- 

+ 

- 

3 

- 

+ 

+ 

4 

+ 

- 

- 

5 

+ 

- 

X 

6 

+ 

+ 

- 

7 

+ 

+ 

+ 

8 

Volumes  intersecting  planes 

— 

9 
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Table  III  and  the  following  definitions  of  the  parameters  IN,  JN  and  KN 


IN  =  0,  x  <  0 
=4,  x  >  0 
JN  =  0,  y  <  0 
=2,  y  >  0 
KN  =  0,  z  <  0 

=  1,  z  >  0  (44) 

are  used  as  an  algorithm  for  computation  of  octants  with  coordinate  compo¬ 
nents.  The  octant  of  a  point  is  defined  as  9  if  either  x,  y,  or  z  =  0. 

The  following  outline  delineates  the  steps  necessary  to  determine  the 
octant  numbers  of  the  elemental  volumes. 

Table  III 
Octant  Logic 


Octant  No. 

IN 

JN 

KN 

Sum  =  IN+JN+KN+1 

1 

0 

0 

0 

1 

2 

0 

0 

1 

2 

3 

0 

2 

0 

3 

4 

0 

2 

1 

4 

5 

4 

0 

0 

5 

6 

4 

0 

1 

6 

7 

4 

2 

0 

7 

8 

4 

2 

1 

8 

Sphere 


The  octant  of  the  center  of  the  sphere  is  determined;  then,  the  following 
differences  are  checked  for  negative  values: 

,  -D  (1)  I  -  R 
1/D  a)  I  -R 

I  ZD  (!)  I  -  R  (45) 

If  there  are  no  negative  numbers,  the  octant  number  of  the  center  point  is 
used.  If  there  are  negative  numbers  octant  9  is  used. 

Hemisphere 


The  octant  numbers  of  the  two  points  which  define  the  hemisphere  are 
determined.  If  the  two  octant  numbers  are  identical,  differences  in 
Equation  45  are  checked  for  negative  values.  If  there  are  no  negative 
numbers,  the  octant  number  of  the  two  points  is  used.  Otherwise,  octant  9 
is  used. 

Cylinder 


The  octant  numbers  of  the  two  points  on  the  center  line  of  the  cylinder 
are  determined.  If  the  two  points  have  the  same  octant  numbers,  the  differ¬ 
ences  in  Equation  45  are  checked  for  negative  values  for  points  1  and  2 
(Figure  4).  If  there  are  no  negative  numbers,  the  octant  numbers  of  the  two 
points  are  used.  Otherwise,  octant  9  is  used. 

Cone 


The  octant  numbers  of  the  points  (1,  2)  for  the  cone  c.nd  (1,  3)  for  the 
truncated  cone  are  determined.  If  the  two  points  have  the  same  octant,  the 
following  differences  are  checked: 

|  xD  (I)  |  -  R (I) 

I  yD  (I)  I  -  R  (i) 

I  ZD  (D  I  -  R(I)  (46) 
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where,  for  the  regular  cone 

1  =  1  and  R(l)  =  radius  of  the  base  and  for  the  tru  .  ated  cone 

and  for  the  truncated  cone 

I  =  1,2  with  R  ( 1 )  =  radius  of  the  base 
R(2)  =  truncation  radius 

If  there  are  no  negative  differences,  the  octant  number  of  point  1  is  used. 
Otherwise,  octant  9  is  used. 

Ellipsoid 


The  octant  of  point  4  of  center  of  the  ellipsoid  is  determined.  The 
following  differences  are  checked  for  negative  values: 


I  *D  (4)  j  -  R 

I  yD  (4)  j  -  R 
lZD  <4)i-Rm 


m 


m 


(47) 


where 


R  =  maximum  major  axis 
m 

If  all  differences  are  positive,  the  octant  of  point  4  is  used.  Otherwise, 
octant  9  is  used. 


Hexahedron 


The  octants  of  the  eight  points  defining  the  hexahedron  are  determined. 
If  they  are  identical,  the  computed  octant  is  used;  otherwise,  octant  9  is 
used. 


This  concludes  the  analysis  required  to  determine  the  individual  shield 
volume  path -length  computations.  The  computed  path  lengths  are  now 
converted  from  inch  units  to  areal-density  units. 
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7.  AREAL-DENSITY  FUNCTIONS 

The  previous  subsections  discussedtheprocedure  for  selecting  the  source 
particle  directions  and  penetration  path  lengths  through  the  structural 
materials.  Now  the  incident  source  particles  must  be  traced  from  the 
shield  external  surface  to  the  dosimeter.  Also,  the  path  lengths  are  converted 
from  length  to  areal-density  units. 

Two  forms  of  areal  density  are  computed  for:  (1)  electron  primary  and 
(2)  heavy  charged  primary  dose  computations. 

a.  Electron  Areal-Density  Function 


For  the  electron  computations,  the  volume  path  lengths  (Pij)  are 
ordered  according  to  decreasing  distances  to  the  volumes  (Qij)  and  trans¬ 
formed  to  an  areal-density  array  (A^g).  Therefore 


P.  P.. 
J  U 


P. 

xr 


A.  »  p  P. 
in  m  lm 


(48) 


where 


Q.. 


>  Q. 


ir 


and 


>Q. 


lm 


P„  =  path  length  for  the  ith  ray  and  the  jth  material 

Q„  =  distance  to  the  jth  path  length  of  the  ith  ray 

p.  =  material  density  of  the  jth  path  length  for  the  ith  ray 
J 

A i^  =  path-length  areal  density  in  (grams /cm2)  units  for  the  ith  ray 
and  the  £th  path  length 
n  =  number  of  path  lengths  for  the  ith  ray 

b.  Heavy  Charged  Particle  Areal-Density  Function 


It  is  desirable  to  convert  the  path  lengths  for  the  various 
equivalent  path  lengths  in  a  standard  material,  which  will  have 


materials  to 
the  same 
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effective  heavy  charged  particle  ionization  transport  properties.  The  equiv¬ 
alent  path  lengths  can  be  computed  without  ordering  the  path  length  according 
to  distance. 

The  conversion  to  a  standard  material  permits  propagation  of  the  total 
spectrum  from  incidence  to  emergence.  This  eliminates  the  necessity  of 
considering  individual  spectral  energies  and  their  propagation  through  each 
material  encountered. 

If  the  range  (R)  versus  energy  (E)  is  approximated  by  the  empirical 
equation 

R  =  SE*1  (49) 

the  energy  transport  equation  is 

R(E  )  =  R(E  )  +  x  (50) 

o  1 

where 

Eq  =  charged  particle  incident  energy 
E^  5  charged  particle  emergent  energy 
x  =  shield  areal  density  (grams /cm2) 

From  Equations  49  and  50 

‘.O'  =  *sEin"+*.  (51) 


where  subscriot  s  refers  to  a  standard  material.  Also 


^A  ^A 

6  A  Eo  6A  E  +  XA 


(52) 


and  A  refers  to  an  actual  material.  Therefore 


xs  =  x.  ,E0-  El> 


XA  ~  XA  {E0'  El^ 


(53) 


E^  is  eliminated;  then 


XS  =  X8  (XA-  V 


(54) 
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However,  to  the  first  order,  xg  can  be  determined  as  a  function  of  as 
follows : 


For  Ex  =  0 


x  =  6  E. 

s  s  0 


XA  =  6AE0 


xa\a 


So,  to  the  first  order 


x  =  5 
s  s 


ft) 


The  error  in  this  approximation  increases  as  a  function  of  increasing 
emergent  energy  and  with  smaller  atomic  weights.  However,  it  is  a  good 
approximation  for  most  practical  applications* 

With  systematic  or  random  selection  of  the  proton  incident  direction, 
an  accumulative  distribution  function  of  the  solid  (Sa)  as  a  function  of  a  real 
density  (T)  is  computed,  where 

Sa  (T)  5  fraction  of  the  total  unit  sphere  with  areal  densities 
less  than  T 

If  Sa  (T)  is  determined  as  a  polynomial  least-square  fit,  then 


S  s  S  C.  Tl 
a  jLj  i 
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Thr  fraction  of  solid  angle  (Sa)  less  than  an  arbitrary  thickness  (T) 
can  be  utilized  for  future  dose  calculations  with  the  same  material  and 
stiuctural  configuration  and  dosimeter  location.  The  areal-density  range  of 
values  may,  for  example,  be  0.  1  to  40  grams  per  centimeter  square.  How¬ 
ever,,  the  smaller  thickness  range  would  contribute  the  larger  portion  of  the 
dose  or  emergent  flux.  Therefore  it  is  necessary  to  place  larger  weight  or 
significance  on  the  smaller  thicknesses.  This  is  done  by  using  more  joints 
for  the  least-square  fit  of  Sa(T)  at  smaller  values  of  T. 

For  discrete  selection  of  T  values  (see  Figure  17)  let  the  value  of  the 
weight  function  for  nth  thickness  be 


R 

n 


n6R 


R  .  +  (1-R  .  ) 
min  mm 


then 


T 

n 


T  .  + 

min 


T  -  T  . 
max  min 


IT 


-1 

cos 


(1-U  ) 
n 


where 


u 

n 


/  h&R  -  R  .  \ 
I  _ min  j 

\  1  -  R  .  / 

'  min  ' 


T 


The  range  of  n  values  is 


• 

/R  .  > 

/, 

n\ 

IPO 

B 

[  min 

\  6R  i 

)  +1 

to 

Ipofci 

r)j 

where 

IPO  (X)  =  Integral  part  of  (X) 


(60) 


(61) 


(62) 


In  general,  the  weight  function  shape  parameter  (t  )  should  be  selected 
inversely  proportional  to  the  minimum  areal  density  of  the  distribution 
function.  Section  III  describes  the  computer  program  and  correlates  its 
computational  operations  with  the  analysis. 
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SECTION  III 


MODIFIED  ELEMENTAL  VOLUME  DOSE  PROGRAM 


A  functional  block  diagram(Figure  18)  depicts  the  operational  functions 
of  the  ray  tracing  subroutines  of  the  modified  EVDP  (MEVDP).  These 
operations  are  correlated  with  the  subroutines  in  Table  IV.  This  table 
illustrates  the  functional  parts  of  the  MEVDP.  The  program  is  organized 
into  three  major  functional  parts:  (1)  data  input  and  transformation,  (2)  ray 
tracing,  and  (3)  areal-density  computations.  The  data  input  and  transforma¬ 
tion  section  of  the  program  reads  the  astronaut  and  spacecraft  input  data  in 
their  respective  coordinate  systems.  Then  the  astronaut  data  are  modified 
to  account  for  limb  rotation.  Subsequently,  both  sets  of  data  are  transformed 
to  a  coordinate  system  positioned  at  the  dosimeter  location.  The  trans¬ 
formed  elemental  volume  data  are  modified  by  matrix  or  coordinate  axes 
rotations.  The  ray  tracing  section  of  the  program  computes  the  ray  direc¬ 
tion  cosines  with  two  options:  (1)  random  and  (2)  systematic  selection.  The 
random  option  is  desirable  for  shielding  configurations  with  a  relatively 
large  degree  of  spatial  inhomogeneity.  However,  the  systematic  selection 
option  is  suitable  for  most  applications.  A  North  American  Rockwell  com¬ 
parison  of  these  options  is  published  in  Reference  1.  The  third  part  of  the 
program  uses  the  data,  which  will  derive  from  the  other  two  parts  to  calcu¬ 
late  the  desired  areal-density  data  be  used  in  AFWL  secondary  and  primary 
nuclear  radiation  programs. 

Further  explanation  is  required  to  amplify  and  reiterate  the  purpose 
of  the  coordinate  systems  introduced  in  Figure  18. 

The  spacecraft  geometrical  structural  configuration  is  originally 
defined  by  elemental  volume  coordinates  plus  parameters  in  the  ABCS  (Fig¬ 
ure  18).  In  order  to  substantially  decreas*  the  amount  of  ray  tracing  or 
tracking  required,  the  elemental  volumes  are  classified  according  to  their 
octant  location  in  a  coordinate  system.  However,  if  the  volume  octants  are 
selected  relative  to  the  ABCS,  there  will  be  unnecessary  tracking.  This 
results  from  the  fact  that,  for  isotropic  source  ray  tracing  calculations,  the 
rays  are  directed  toward  the  dosimeter  which  is  generally  displaced  from 
the  origin  of  the  ABCS.  Therefore,  the  ray  could  traverse  four  octants. 
However,  if  the  geometry  is  transformed  to  a  coordinate  system  centered  at 
the  detector  (DSCS),  the  rays  traverse  only  one  octant.  Tracks  along  the 
axes  are  considered  as  special  cases.  The  axes  of  the  ABCS  and  the  DSCS 
have  the  same  vectorial  direction  and  therefore  differ  only  by  a  vectorial 
displacement. 


For  several  elemental  volumes  it  is  convenient  to  have  a  detector 
coordinate  system  with  one  or  more  of  its  axes  directed  along  the  reference 
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Table  IV 

Functional  Description  of  the  MEVDP  Subroutines 


r 


i 

3 
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(Continued) 

if  the  MEDVP  Subroutines 


areal-density  entrance  path  length  and  material  type 


Table  IV  (Concluded) 

Functional  Description  of  the  MED  VP  Subroutines 


Figure  18  is  located  on  page  215. 


Figure  18.  Functional  Diagram  of  the  MEVDP 
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axes  of  the  elemental  volumes.  Since  the  new  detector  system  has  th  ;  same 
origin  as  the  DSCS,  the  ray  direction  cosines  can  be  converted  to  the  new 
system  by  matrix  rotations.  The  new  system  is  designated  the  Rotated 
Detector  Coordinate  System  (RDCS).  In  order  to  select  the  rotation  matrix 
transformations  for  the  elemental  volumes,  it  is  necessary  to  determine 
convenient  coordinate  systems  relative  to  the  elemental  volumes  reference 
axes.  The  most  convenient  axes  are  selected  as  those  which  facilitate 
calculations  of  path  lengths  to  the  elemental  surfaces. 

1.  MAIN  PROGRAM 

The  main  program  controls  the  calling  sequence  of  the  major  routines 
for  data  input  and  transformation,  ray  tracing,  and  areal-density  computa¬ 
tions.  For  the  reader's  convenience,  the  discussion  coordinates  the 
FORTRAN  IV  program  listing  (Appendix  II)  with  the  first  section  of  this 
document.  R  I  is  used  to  define  the  reference  statement  I  which  is  nearest 
to  the  operation  being  discussed.  These  references  are  comment  cards,  in 
the  program  listing,  entitled  "Program  Operation  Reference  No.  I,  "  where 
I  is  the  number  in  R  I  .  Reference  comment  cards  before  DO  loops  refer 
to  the  operations  performed  in  the  indicated  loop.  Scratch  disk  storage 
devices  are  designated  as  SYSDA.  The  main  program  rewinds  tape  units  8 
and  9  and  reads  the  parameters  N10GEN  and  CONTEN.  CONTEN  is  a 
descriptive  title  and  N10GEN  specifies  that  the  input  geometry  is  read  as 
cards  or  from  tape  unit  10.  At  Rl,  subroutine  GENTAP  is  called  to  generate 
a  card  image  of  the  geometry  data  on  tape  unit  10  if  N10GEN  >0.  At  R2, 
the  input  parameter  OPTN  transfers  to  R3  to  input  a  new  set  of  geometry 
or  continues  and  calls  subroutines  CLOCK  and  TRSHLD.  CLOCK  is  used 
to  compute  the  elapse  time  in  TRSHLDtwhich  transforms  the  geometry  data. 
Then  routines  ESDOSE  and  ORDER  are  called  to  (1)  perform  ray  tracing 
and  areal -density  computations,  and  (2)  order  the  areal-density  data  accord¬ 
ing  to  ray  numbers  and  decreasing  distance  from  the  dosimeter  location. 
Subroutine  CLOCK  uses  the  library  function  SECOND  to  compute  current 
time  in  seconds,  and  converts  the  time  to  minutes. 

2.  SUBROUTINE  GENTAP 

The  function  of  this  subroutine  is  to  read  geometry  card  input  data  and 
generate  two  files  containing  (1)  the  basic  spacecraft  geometry  and  (2)  the 
simulated  astronaut  geometry  on  tape  unit  10  in  card  image.  The  end  of  the 
input  data  is  determined  by  testing  after  R4  for  a  blank  card. 

3.  SUBROUTINE  FILE 

Subroutine  FILE  moves  tape  unit  number  NTAPE  to  the  beginning  of  the 
file  numbered  (NEWF+1).  If  tape  number  NTAPE  has  been  used,  the  old 
file  number,  NOLDF,  of  NTAPE  is  used  with  NEWF  to  avoid  rewinding  when 
tape  NAPE  can  be  advanced  to  the  proper  file  (i.  e.  ,  NSKIP  >  0  at  R5). 


4,  SUBROUTINE  TRSHLD 
a.  Geometry  Data  Transformation 

This  subroutine  uses  input  geometry  data,  which  define  the  elemental 
volumes  (Figure  4)  and  the  dosimeter  location.  The  original  coordinates, 
in  the  ABCS,  are  transformed  to  the  DSCS  by  a  translation  (Equation  11); 
then  the  DSCS  and  tracking  ray  direction  cosines  are  transformed  to 
the  RDCS  by  an  equivalent  matrix  rotation  of  the  DSCS,  and  the  octants 
computed  in  the  DSCS  for  each  shield.  Auxiliary  computations  include  the 
ellipsoid  major  axes,  radius  of  the  base  of  the  cone,  etc. 

The  composite  transformed  and  auxiliary  data  for  each  shield  is 
written  on  a  SYSDA  unit.  This  tape  will  be  input  data  for  the  areal-density 
computations.  The  following  discussion  presents  the  definitions,  analysis, 
and  logic  used  to  implement  the  functions  of  TRSHLD. 

The  variables  N  and  NSO  are  set  to  zero,  where 


N  =  a  counter  which  is  equal  to  the  number  of  the  elemental  volume 
being  transformed 

NSO  =  number  of  shields  in  the  Ith  octant  in  the  DSCS 


Then  the  following  input  data  are  read: 


XDET(I)  =  detector  coordinates  in  the  ABCS  and  1=1,  2,  3  refer 
to  the  x^»  YA»  ZA  components 

NS  =  number  of  input  elemental  volumes  or  a  positive  integer 
for  one  or  more  volumes,  and  zero  for  no  spacecraft 
volumes 


PRNT  =  printout  control  parameter  (=0,  no  detail  tracking  print¬ 
out  desired;  *  0,  all  printout  required) 

NOMAN  =  total  number  of  elemental  volumes  comprising  one 
astronaut 


NMEN  =  number  of  astronauts  (5 10) 

NLIMB(I)  =  number  of  composite  shields  which  will  be  rotated  for 
the  Ith  astrouaur  and  is  defined  as  unity  for  the  NASA- 
MSC  model  astronaut 


NOMAN  is  redefined  as  the  total  number  of  elemental  volumes 
comprising  the  crew 
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Tape  unit  12  (SYSDA),  which  will  contain  the  transformed  data,  is 
rewound.  Then  the  elemental  volumes  counter,  N,  is  increased  by  one 
and  tested  to  determine  if  all  shields  have  been  transformed.  The  following 
parameters  are  read  as  input  data  for  the  Nth  shield  from  tape  10: 


SSN  = 
ST  = 
MATCD  = 

RPHI  = 
RHO  = 
JMAX  = 
NCOMP  = 


shield  serial  number 
shield  type 

material  code  number  which  defines  the  shield  material 
composition 

radius  of  cylinder  or  radian  half-angle  of  a  cone 

shield  density  (grams  /cm^) 

type  of  ellipsoid  (degree  of  truncation) 

the  number  of  elemental  volumes  comprising  a  composite 
shield  (£10) 


As  many  as  ten  of  the  elemental  volumes  may  be  combined  to  form 
more  complex  shields,  which  are  designated  as  composite  shields.  The  void 
(negative)  elemental  volume  components  are  used  to  void  the  spatial  region 
they  occupy.  The  composite  shield  permits  imbedding,  in  which  the 
positive  (solid)  shield  may  contain  all  or  a  portion  of  a  void  shield.  Also, 
the  solid  shield  components  may  have  different  material  types  and  densities. 
Portions  of  several  of  the  positive  shield  components  may  contain  portions 
of  a  single  void  shield.  This  versatility  results  in  a  very  significant  improve¬ 
ment  in  the  geometrical  representation  capability  and  reduces  the  magnitude 
of  the  shield  synthesis  effort. 


The  variable  TG  is  used  to  define  tape  unit  10.  The  shield  types  (ST) 
are  defined  as  follows: 


Type  of  Shield 

Hollow  Shield 

Solid  Shield 

Hexahedron 

0 

1 

Cylinder 

2 

3 

Sphere 

4 

5 

Hemisphere 

6 

7 

Cone 

8 

9 

Truncated  Cone 

10 

11 

Ellipsoid 

12 

13 
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JMAX  is  defined  as  follows,  and  is  applicable  only  to  elliptical  volumes: 


JMAX 

Type  of  Ellipsoid  Truncation 

4 

Non-truncated 

+5 

Single  truncation  with  lower  portion 

Eliminated  with  respect  to  z^  axis 

-5 

Single  truncation  with  upper  portion 

Eliminated  with  respect  to  +zr  axis 

6 

Double  truncation 

The  next  portion  of  the  program  through  R6  uses  ST  to  determine  the 
number  of  coordinate  points  used  to  define  the  shield.  From  Figure  4,  the 
number  of  coordinates  as  a  function  of  shield  types  is  as  follows: 


ST 

Number  of  Shield  Coordinates 

0,  1 

8 

2  through  9 

2 

10,  11 

3 

12.  13;  JMAX  =  4 

4 

=  ±4 

5 

=  6 

6 

After  the  program  logic  is  used  to  compute  the  number  of  coordinates, 
the  coordinates  [X(I,  J)]  are  read  as  input  data  from  tape  unit  TG.  I  and  J 
refer  respectively  to  the  point  number  and  coordinate  components  (x^,  yA» 
zA).  NP  is  equivalent  to  the  maximum  number  of  coordinates  for  a  given 
shield  type. 
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Between  R7  and  R8,  the  shield  coordinates  are  converted  from  the 
two-dimensional  array  X(I,  J)  to  a  one -dimensional  array,  XTRAN(IJ),  with 
the  following  equivalence: 


I 

J 

IJ 

1 - NP 

1 

1 - -  NP 

1 — -NP 

2 

NP+1 - -2NP 

1  — ►NP 

3 

2NP+1 — -3NP 

Since  ST  =  0  is  not  permitted  in  the  subsequent  GO  TO  logical  statement,  it 
will  be  incremented  by  1  and  defined  as  1ST.  The  logical  statement  after 
R9  transfers  control  to  subroutines  which  perform  the  transformations 
for  the  required  shield  type  (see  Functional  Block  8  of  Table  IV). 

The  translation  and  rotation  transformations  are  performed  for  the 
Nth  shield  as  follows: 

Translation 

Since  the  translation  to  the  OSCS  is  performed  in  a  similar  manner 
for  each  shield,  it  will  be  shown  only  for  the  hexahedron.  Through  equiva¬ 
lence  XTRAN  is  redefined  as  XHEX.  From  Equation  11,  the  translation 
from  the  ABCS  to  the  DSCS  is 

XHEX  (I,  J)  =  XHEX(I,  J)  -  XDET(J) 

Another  example  of  this  translation  is  RIO  of  subroutine  CYL.NDR. 

Rotation 


For  the  hexahedron,  the  rotation  matrix  is  based  on  the  vectorial 
direction  of  the  tracking  ray,  which  is  not  yet  defined-  Therefore,  no  rota¬ 
tion  is  done  for  this  shield  in  the  geometry  data  transformation  subroutine. 

For  the  other  elemental  volumes,  the  generalized  matrix  of  rotation 
(Equation  34)  is  redefined  by  Mxi(8/)  •  Mza(^  )  =  A(K,  J),  K  =  1,3  and 
J  =  1,3.  — 


The  trigonometric  functions  of  the  angles  of  the  rotation  matrix  are 
computed  as  follows:  a  vector  distance  between  points  2  and  1  is  defined  by 
the  components 

D X(K)  =  X(2,K)  -  X(1,K);  K  =  1,3 
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From  Eqs.  35  and  36 
TANPHI  =  tan  4> 

COSPHI  =  cos  0 
SINPHI  =  sin  0 

The  magnitude  of  the  vector  is  CL  and 
COSTH  =  cos  6 
SINTH  =  sin6 

To  account  for  the  signs  of  the  trigonometric  functions  as  function  of  quadrant, 
the  logical  equivalents  to  Equations  35  and  36  are  executed  between  Rll  and 
R12  for  the  cylinder  subroutine  (CYLNDR).  The  same  logic  is  also  used  for 
the  hemisphere  and  cones. 

Now  the  rotation  matrix  is  defined  and  computed  from  Equation  34 

where  XR  =  X  (RDCS),  XD  =  X(DSCS) 

with  the  equivalent  XR  =  A  -XD  >or 

3 

XR(I,  K)  =  ^  A(K,  J)XD(I,  J) 

J=1 


where  1=1,  NP  points  and  K  refers  to  the  3  components 

Rotation  of  the  ellipsoid  from  the  DSCS  to  the  ROCS  requires  a  special 
matrix  (see  Subroutine  ELIPSD).  The  Rk  vectors  defining  the  major  axes 
are  computed  as  follows: 


The  magnitudes  and  direction  cosines  of  the  vectors  are  DIST(K), 
where  K  refers  to  the  three  axial  vectors  in  the  DSCS  and  R13  of 
ELIPSD,  and 


DIR(K,  J) 


XD(K,  J)  -  XD(4,  J) 
DIST(K) 


J  =  1,  2,  3  and  K 


1,  2.  3 
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From  Equation  31,  the  FORTRAN  equivalent  for  the  rotation  matrix  is 
DIR  (I,  J),  and  the  coordinates  in  the  RDCS  are 


3 

XR(I.K)  =  2  DIR(K,  J)  •  XD(I,  J)  where  *(I)DJ  =  DIR(I,J) 

J=1 


The  lengths  of  the  ellipsoid  axes  (a,  b,  c), which  are  defined  as  SMMA(K),are 
computed  in  the  rotated  system  as 

SMMA(K)  =  XR(K,  K)  -  XR(4  K),  K  =  1,  3  axes 

Since  the  standard  astronaut  is  defined  in  a  unique  coordinate  system, 
special  transformations  are  required. 

The  standard  man  is  a  NASA-MSC  model  which  is  used  for  shielding 
calculations  (Reference  2).  Figure  19  shows  a  diagram  of  the  man  with  his 
cylindrical  limbs  designated  Cl  through  C6.  The  Standard  Man  Coordinate 
System  (SMCS)  is  a  right-handed  orthogonal  set  with  the  x  axis  directed 
from  the  back  toward  the  front  of  the  standard  man.  The  man's  feet  are 
represented  by  the  hexahedrons  H6  and  H8.  The  correspondence  between 
the  NASA-MSC  man's  limbs  and  the  North  American  Rockwell  Space  Division 
(NR-SD)  shield  serial  numbers  (SSN)  is  as  follows: 


NASA-MSC 

Man's  Limb 

NR-SD 

Shield  Serial  Number 

Cl 

9010 

C2 

9011 

C3 

9012 

C4 

9013 

C5 

9014 

C6 

9015 

C7 

9016 

C8 

9017 

H6 

9018 

H8 

9019 
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Figure  19.  NASA-MSC 
Standard  Man  Model 
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In  Figure  19  the  standard  man  is  in  the  standing  position.  Therefore, 
the  limbs  and  feet  must  be  rotated  in  order  to  place  him  in  reclining  or  sitting 
positions  relative  to  the  SMCS.  Additional  rotations  and  translations  will  be 
required  to  transform  the  man's  numerical  shield  coordinates  to  the  ABCS. 
The  man's  limbs  are  rotated  with  a  vector  matrix.  Let  VR  and  Vp  be  the 
vector  position  of  a  point  (R)  to  be  rotated  relative  to  a  fixed  point  (F).  Then 
(Vr-Tf)  is  the  vector,  representing  the  limb,  which  is  to  be  rotated.  The 
new  vectorial  position  of  point  R  is  VRj  where 

7ri  =  <VVi+7f  (64> 

If  (VR  -Vp )  is  rotated  through  an  angle  e  about  the  y  axis,  using  an  equation 
developed  by  Rodriguez,  where 


E,  =  V-V_ 
1  R  F 

E  =  T 


ix+jy+kz 


p  =  e 


E  =  (V  -V  ) 

2  '  R1  F'l 


E2  =  Ej  +  (cos  p-1)[ei  -  (E-Ej)  E  J  +  sinp  E  xEj  (65) 

In  general  the  vector  Ej  is  rotated  about  E  through  an  angle  p  in  the  direction 
of  a  right-hand  screw  directed  along  E.  If  the  original  angle  between  E  and 
Ui  is  constant  during  the  rotation,  the  vector  Ej  is  transformed  to  E2. 


Since 


E  *  Ej  =  y 
ExEj  =  iz-kx, 


in  matrix  notation 


<VR-Vl 


cos  c  0  sin  c 


-sin «  0  cos  t 


(VR-  V 


-  a  <VR-  V 
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Equations  65  and  66  are  used  to  compute  the  shield  coordinates  Vpi  after 
the  limb  and  feet  rotations.  For  the  cylinders  (Figure  19).  point  1  is  defined 
as  the  point  with  the  larger  z  value.  The  fixed  points  for  limb  rotations  will 
always  be  point  one  of  a  cylinder.  The  sequence  of  rotations  is  defined  in 
Table  V  in  terms  of  the  fixed  point,  rotated  point  and  angle  of  rotation.  The 
limb  rotated  man  is  transformed  to  the  ABCS  by  a  combination  of  translation 
and  rotation.  In  Figure  20,  points  1  through  4  are  used  to  define  the  desired 
position  and  orientation  of  an  astronaut  relative  to  the  ABCS. 

If  point  P  of  one  of  the  man's  elemental  volume  components  is  defined 
as  %  relative  to  the  SMCS,  then  its  coordinates  can  be  defined  relative  to 
the  ABCS  as  Ra  where 


(67) 


Equation  67  represents  the  required  translation.  However,  to  effect  the 
translation  Rm  must  be  transformed  to  the  ABCS.  Let 


(68) 


where  a  refers  to  summation  notation,  defined  after  Equation  27,  and  x  and 
y  are  coordinate  components  in  the  SMCS  and  ABCS  respectively.  From 
Figure  20  _ 

T  =  lA  ;  or  =  1.  3  ;  j  =  1,  3  (69) 

Mj  Jo  \ 

where 


(70) 


and 


j  =  coordinate  axes  points  (see  Figure  20)  where  (1,  2,  3,  4)  corre¬ 
spond  to  the  (x,  y,  z  axes  and  the  origin) 

m  refers  to  a  point's  coordinate  components 

For  the  next  derivation  the  summation  convention  is  further  defined  as 
follows;  if  in  some  expression  a  certain  index  occurs  twice,  it  shall  mean 
that  this  expression  is  summed  with  respect  to  that  index  for  all  admissible 
values  of  that  index. 
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Table  V 

Fixed  and  Rotated  Points  and  Angles  for  Astronaut  Limb  Rotations 


Rotation 

Number 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12-19 

20 

21-28 

29 

30 

31 
32-39 

40 

41-48 


Fixed  Point 


CCN 


CPN 


CCN 


Rotated  Point 


CPN 


Angle 
of  Rotation 


1 

1 

1 

2 

3 

3 

3 

4 

5 
5 
5 

5 

6 
6 
7 
7 
7 

7 

8 
8 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

2 

2 

2 

3 

4 
4 

4 

5 

6 
6 
9 
6 
9 

7 

8 
8 

10 

8 

10 


2 

1 

2 

2 

2 

1 

2 

2 

2 

1 

2 

1-8 

2 

1-8 

2 

1 

2 

• -8 
2 

1-8 


1 

1 

1 

2 

3 

3 

3 

4 

5 
5 
5 

5 

6 
6 
7 
7 
7 

7 

8 
8 


#CCN  a  Composite  shield  component  number  (SSN  9010  —  CCN=1  etc.  ) 
CPN  *  Component  shield  point  number 
Angle  I  ■  «  j 


Summarizing,  Equations  67,  70,  and  72  will  transform  the  man's  elemental 
volume  coordinates  to  the  ABCS.  The  transformation  of  the  crew  shields 
to  the  ABCS  are  executed  in  subroutines  FANTOM  R14  and  LIMROT  (limb 
rotations),  after  which  (R15  of  TRSHLD)  the  astronaut  elemental  volumes 
are  transformed  to  the  RDCS. 

b.  Octant  Determination 

Hexahedron 

The  octant  of  the  position  of  the  8  apexes  is  computed  in  the  DSCS. 

If  any  component  of  the  coordinates  of  any  point  is  zero,  the  octant  is  defined 
as  9-  That  is,  the  figure  may  be  in  more  than  one  octant. 

Equation  44  and  Tables  2  and  3  are  implemented  with  the  following 
equivalent  variables: 

XHEX(I,  J)  =  (X  ,  Yt,  Z  );  J  = 

1  I  = 

IN(J,  K)  =  (INj,  JNJf  KNj);  K  = 

J  — 

Then  the  octant  of  the  Jth  point  is  IN  (J,  1)  =  IN(J,  1)  +  IN(J,  2) 

+  IN(J ,  3)  +  1.  If  any  IN(J,  1),  J  =  1,  8  are  different,  then  the  octant 
number  is  9.  Otherwise  the  octant  is  IN(J,  1). 

The  other  elemental  volumes,  or  volumes  which  contain  them,  are 
checked  for  intersections  with  the  coordinate  planes  in  the  DSCS. 

A  slight  variation  of  Equations  45  to  47  is  used  to  determine  if  each 
volume  intersects  any  coordinate  planes.  The  method  is  as  follows: 

Equation  45  is  replaced  by 

XD  (1)  ±R 

yd  (1)  ±R 

ZD  (1)  ±R 

If  signs  of  [xD  (1 )  +  r]  and  [xD(l)  -  R  ]  and  the  similar  terms  of  the 
other  components  are  the  same,  there  is  no  intersection  with  the  coordinate 
planes. 


1.2.3 

1,2,.. ..,8 

1.3 

1,  2,  ...  ,8 
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This  method  is  incorporated  in  the  program  as  follows: 

The  signs  of  |x£>  (1)  +  r],  [xq  (1)  -  R  j  are  determined  with  the  aid  of 
FORTRAN  IV  library  function  SIGN.  SIGN  is  used  to  define  Sj  and  S2 
as  follows: 


Sx  =  SIGN 


S„  =  SIGN 


Xd(I,K),  |xd(I,K)  +  r] 

|sign  of  Xd(I,K)  +  R  ]*  j  Xd(I,K)|,  and 
Xd(I,K).  [xd(i,k)-r] 


If  Sj  +  S2  =  0  for  either  coordinate  component,  the  elemental  volume 
intersects  a  coordinate  plane.  Otherwise  the  octant  of  the  coordinate  of 
the  center  of  volume  is  used.  For  the  ellipsoid  R  is  replaced  by  the  largest 
major  axis. 


Auxiliary  Computations 

The  radius  of  the  sphere  and  hemisphere  is  computed  as  the  magnitude 
of  the  vector  between  points  2  and  1  of  the  volume. 


The  radius  (RCON)  of  the  base  of  the  cone  is  computed  as  follows: 

PHICON  =  RPHI  is  defined  as  the  radian  half-angle  of  the  cone 

The  magnitude  (CL)  of  the  vector  between  points  2  and  1  is  used  to 
compute  the  radius  by 


RCON  =  CL  tan  (PHICON)  —  CL  *  TANCON 


The  radii  of  the  two  bases  of  the  truncated  cone  are  computed  as  follows 
(Figure  8): 

The  distance  (DZ)  from  point  2  to  3  is  determined  with  the  RDCS 
points  by 


DZ  =  XR(2,  3)  -  XR  (3,  3) 

then  the  radii  of  the  bases  at  points  3  and  1  (RT,  RB)  are 

RT  =  DZ  *  TANCON 
RB  =  RCON 
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This  concludes  the  algebraic  and  logical  operations  necessary  to 
perform  the  data  transformations,  octant  computations,  and  auxiliary 
computations.  All  information  for  each  shield  is  stored  in  a  single  sub¬ 
scripted  array  called  SHLDD  at  R16  in  TRSHLD.  The  maximum  number  of 
locations  used  for  a  shield  is  38,  i.  e.  ,  for  the  ellipsoid.  Although,  for 
convenience,  all  of  the  SHLDD  arrays  contain  38  words,  some  words  are 
undefined.  NSO(I)  is  used  as  a  counter  to  determine  the  number  of  shields 
in  the  Ith  octant.  The  components  of  the  SHLDD  array  are  as  follows: 

In  general  the  array  contains: 

SHLDD(l)  =  OCT  s  octant  in  which  shield  lies 

SHLDD(2)  =  ST  s  shield  type 

SHLDD(3)  =  SSN  *  shield  serial  number 

SHLDD(4)  =  MATCD  m  material  code 

SHLDD(5)  =  RHO  2  density  of  materials  (gm/cm.  3) 

SHLDD(6)  =  NP  *  no.  sets  of  coordinates 

SHLDD(7)  =  RPHI  *  radius  (inches)  or  jhalf-angle  of  cone  (radians )J, 
whichever  applicable 

SHLDD(8  -  37)  =  XR(I,  J),  1*1,2 .  NP  and  J  =  1,  2,  3,  = 

coordinates  (inches)  in  rotated  detector  system 
and  other  information 

SHLDD(38)  =  NCOMP  =  number  of  composite  shield  components 
FOR  SPECIFIC  SHIELD  TYPES: 

1.  SHLDD(2)  =  0  or  1  for  hollow  or  solid  hexahedron,  respectively 

SHLDD(6)  =  8  =  no.  corners  of  hexahedron 

SHLDD(8  -  15)  =  XHEX(1, 1);  SHLDD(16-23)  =  XHEX(1,2): 

I  =  1,2,..  .,8 

SHLDD(24-31)  =  XHEX(I,  3);  I  =  1 , 2, .  .  . ,  8 

SHLDD(32-37)  contain  zeros,  not  needed  for  the  hexahedron 

NOTE:  Locations  not  mentioned  contain  information  described 
in  general  for  all  shields. 


65 


2.  SHLDD(2)  =  2  or  3  for  cylinder  (hollow  or  solid) 

SHLDD(6)  =  2  sets  of  coordinates  for  the  centers  of  the  bases 
SHLDD(7)  =  RPHI  =  radius 

SHLDD(8-9)  =  XR(I,  1);  SHLDD(lO-ll)  =  X(1, 2);  SHLDD(12-13) 

=  X(I,  3);  1=  1,2 

SHLDD{14-22)  =  Rotation  matrix  A(3,  3)  where  the  order  is: 

A(1 ,  1 ),  A(2,  1 ),  A(3,  1),  A(l,  2),  A(2,2),  A(3,2),  A(l,  3), 

A (2,  3),  A(3,  3) 

3.  SHLDD(2)  =  4  or  5  for  sphere  (hollow  or  solid) 

SHLDD(6)  =  2  sets  of  coordinates  where  the  first  is  the  center 
of  the  sphere  and  the  second  is  a  point  on  the  surface 

SHLDD(7)  =  Radius 

SHLDD(8-9)  =  XD(I,  1);  SHLDD(lO-ll)  =  XD(I,2); 

SHLDD(12 -1 3)  =  XD(I,  3),  1=1,2 

SHLDD(14-37)  contain  zeros 

4.  SHLDD(2)  =  6  or  7  for  hemisphere  (hollow  or  solid) 

SHLDD(6)  =  2  sets  of  coordinates  where  the  first  is  the  center 
of  the  base  and  the  second  is  on  the  z  axis  of  the  hemisphere 
and  on  the  surface 

SHLDD(7)  =  Radius 

SHLDD(8-9)  =  XR(I,  1);  SHLDD(lO-ll)  =  XR(I,2); 

SHLDD(12  -1 3)  =  XR(I,  3),  1=1,2 

SHLDD(14-22)  =  Rotation  matrix  A;  for  order  of  storage  see 
cylinder 

SHLDD(23-37)  contain  zeros 
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5. 


SHLDD(2)  =  8  or  9  for  cone  {hollow  or  solid) 

SHLDD(6)  =  2  sets  of  coordinates  where  the  first  is  at  the  center 
of  the  base  and  the  second  is  at  the  apex 

SHLDD(7)  =  PHI  =  half -angle  at  apex  in  radians 

SHLDD(8-9)  =  XR(1, 1);  SHLDD(lO-ll)  =  XR(I,2); 

SHLDD(12-13)  =  XR{I,  3),  1=1,2 

SHLDD(14-22)  =  Rotation  matrix  A  -  order:  same  as  cylinder 
SHLDD(23)  =  Radius  of  base 
SHLDD(24-37)  contain  zeros 

6.  SHLDD(2)  =  10  or  11  for  truncated  cone  {hollow  or  solid) 

SHLDD(6)  =  3  sets  of  coordinates  where  the  first  is  at  the  center 
of  the  base,  the  second  is  at  the  imaginary  apex  and  the  third 
is  on  ;he  axis  of  the  cone  at  the  plane  of  truncation 

SHLDD(7)  =  PHI  =  half-angle  of  cone  at  imaginary  apex  in 
radians 

SHLDD(8-10)  =  XR(I,  1);  SHLDD(11-13)  =  XR(I,2); 

SHLDD(14-16)  =  XR(I,  3),  1=1,2,  3 

SHLDD(17-25)  =  rotation  matrix  A  -  order  of  storage  same  as 
for  cylinder 

SHLDD(26)  =  Radius  of  base 

SHLDD(27)  =  Radius  of  cone  at  plane  of  truncation. 

SHLDD(28-37)  contain  zeros 

7.  SHLDD(2)  =  12  or  13  for  ellipsoid  (hollow  or  solid)  of  0,  1,  or  2 

truncations 

SHLDD(6)  =  JMAX 
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SHLDD(8-13)  =  XR(I,  1);  SHLDD(14-19)  =  XR(I,2); 

SHLDD(20-25)  =  XR(I,3),  1=  1,  2,  ....  6 

(If|JMAX|<  6,  zeros  are  stored  in  the  unused  locations.  ) 

SHLDD(26-34)  =  direction  cosines  of  ellipsoid  axes,  and  the 
order  of  storage  is  the  same  as  that  of  the  rotation  matrix 
of  the  cylinder 

SHLDD(35 -37)  =  SMMA(K),  K  =  1,  3  =  length  of  the  ellipsoid 
axes  (a,  b,  c) 

At  R17  in  TRSHLD  after  the  hexahedron  shield  data  have  been  trans¬ 
formed,  the  parameter  NX  is  defined.  NX  is  the  number  of  numerical 
coordinates  and  dimension  parameters  as  a  function  of  shield  type.  At  R18, 
the  transformed  dimension  data,  i.  e.  ,  XYZ(I),  are  transferred  to  the  SHLDD 
array.  At  R19,  the  data  for  a  composite  shield  component  are  stored  in  the 
two-dimensional  array  SPN(I,  J),  where  J  refers  to  the  shield  component 
number.  Then,  the  complete  composite  shield  data  are  written  on  SYSDA 
unit  12.  If  all  the  spacecraft  shields  (NOS)  have  been  transformed,  the  end 
of  file  test  after  R20  transfers  to  R14  for  transformation  of  the  crew  shields. 
However,  before  continuing  with  the  program  description,  the  standard  man 
and  the  transformations  required  to  define  the  astronaut  in  the  ABCS  will 
be  delineated. 

5.  SUBROUTINE  FANTOM 

This  subroutine  simulates  up  to  ten  astronauts  in  the  vicinity  of  the 
spacecraft.  The  process  is  accomplished  by  describing  one  standard  man 
in  the  basic  Standard  Man  Coordinate  System  (SMCS)  on  the  input  shield 
geometry  tape  (unit  10).  All  limbs  are  defined  as  one  composite  shield  of 
10  components  and  each  limb  may  be  rotated  through  a  physically  reasonable 
angle  about  the  standard  man's  y  axis. 

Input  arguments  to  this  subroutine  are  as  follows: 

NOMAN.NMEN,  NLB(I)=NLIMB(I)  for  I  =  1,  NMEN,  and  NGT,  where 

NGT  is  equivalent  to  TG  or  tape  number  10 

The  SYSDA  unit  13,  which  will  contain  the  astronaut  geometry  in  the 
ABCS,  is  rewound.  Tape  NGT  is  positioned  at  the  beginning  of  the  second 
file  by  use  of  the  FILE  subroutine  for  each  astronaut.  The  major  loop 
transforms  the  M  =  1,  NMEN  astronauts  to  the  ABCS.  The  number  of  com¬ 
posite  shield  volumes,  which  will  be  rotated  in  the  Mth  man  is  redefined  as 

NLIMB  =  NLB(M) 
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The  direction  cosines  and  coordinate  components  of  Equation  70  are 
defined  in  the  program  as  follows: 

x(j)  =  W(J,M) 
m 

*U>4  =  R4A(J) 

g..  =  L(J,M) 

After  W  and  R4A  are  input,  the  L  matrix  is  computed.  Then  all  required 
rotations  for  the  Mth  man  are  performed.  The  title  record  is  read  from 
tape  NGT  and  reread  because  of  input  requirements  for  TRSHLD  and  this 
subroutine. 

The  first  eight  volumes  on  the  second  file  of  tape  NGT  are  cylinders 
defining  the  man's  arms  and  legs;  the  ninth  and  tenth  volumes  are  hexahed¬ 
rons  defining  the  feet.  The  first  cylinder  is  read  and  stored  in  the  array 
XCMAN(I,  J,  1),  where  j=l,  3  coordinate  components  and  1=1,  2  points  for  the 
cylinder. 

The  number  of  components  (elemental  volumes)  for  rotation  is  defined 
as  NCOMP.  The  10  component  angles  for  rotation  are  read  as  the  array 
EPSLN.  However,  only  eight  of  the  angles  will  be  used  for  the  limb  rota¬ 
tions.  Then  the  header  records  and  coordinates  of  the  remaining  shields  for 
rotation  are  read  as  TITLE(I,  ICP)  and  XCMAN(I,  J,ICP)  where  ICP  denotes 
the  component  number. 

The  subroutine  LIMROT  is  then  called  for  the  limb  rotations.  The 
table  defining  the  fixed  point,  the  point  to  be  rotated,  the  angles  through 
which  the  components  are  to  be  rotated,  and  the  rotation  number  is  used  with 
Equation  65  to  perform  the  limb  rotations.  All  component  coordinates  are 
returned  in  the  array  XRMAN  relative  to  the  SMCS.  Then  XRMAN  is  trans¬ 
formed  by  a  coordinate  rotation  with  Equation  72  and  defined  as  the  variable 
Y.  Finally  a  translation  of  Y  converts  XRMAN  to  the  ABCS.  Now  the 
NCOMP  header  records  and  transformed  coordinates  in  the  ABCS  are  written 
on  SYSDA  unit  13.  This  completes  the  limb  rotations  and  transformation  of 
the  coordinates  (XRMAN)  to  the  vehicle  or  Absolute  Coordinate  System. 

The  other  volumes  of  the  Mth  man,  which  are  not  rotated,  are  trans¬ 
formed  to  the  ABCS.  The  final  coordinates  are  stored  in  the  array  XMT  an4 
the  title  records  along  with  XMT  are  written  on  unit  13  sequentially. 

The  above  process  is  continued  until  M  =  NMEN.  Then  unit  13  is 
rewound,  NGT  is  redefined  as  13  for  use  by  TRSHLD,  and  the  astronaut 
geometry  is  added  to  the  spacecraft  geometry  on  SYSDA  unit  12. 
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6.  SUBROUTINE  LIMROT  (NCOMP,  XCMAN,  XRMAN,  EPSLN) 


This  subroutine  performs  the  rotations  of  the  limbs  of  an  astronaut  in 
the  man  coordinate  system  through  the  required  angles  (EPSLN).  Rotations 
are  made  about  the  y  axis  of  the  standard  man  in  the  SMCS. 

The  subroutine  input  arguments  are  defined  as  follows: 

NCOMP  =  number  of  components  (£10) 

XCMAN  =  unrotated  coordinates  of  all  components 

EPSLN  =  angles  of  rotation 

The  subscripts  of  the  components,  point  numbers,  etc.  ,  are  stored  by  use 
of  the  data  statement  (see  Table  V) 

where 


NCF(I),  I  =  1,  48  =  number  of  the  component  of  the  fixed  point  of  the 

rotated  limb 

NCROT(J),  J  =  1,  48  =  number  of  the  component  of  the  point  which  is 

is  moved  by  limb  rotation 

NCPT(K),  K  =  1,  48  =  number  of  point  which  will  be  moved  by  limb 

rotation 

The  particular  fixed  point  component,  rotation  point  component  and  point 
number  within  the  rotation  component  are  respectively  defined  as 

NF  =  NCF(IR) 

NR  =  NCROT(IR) 

NC  =  NCPT(IR) 

The  IRth  rotation  matrix  o  is  computed  from  Equation  66  and  defined  as 
ALPHA.  Then  ALPHA  is  used  to  transform  the  vector  components  XYZ{K) 
and  the  rotated  coordinates  are  defined  as  XRMAN(NC,  J,  NR).  This 
computation  is  performed  for  the  48  limb  rotations. 

7.  SUBROUTINE  ESDOSE 

This  major  section  of  the  program  computes  the  areal-density  functions. 
One  of  two  methods  will  be  used  for  obtaining  the  direction  cosines  of  a 
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specified  number  of  rays.  One  method  randomly  selects  the  directions  by 
random  sampling  from  an  equally  probable  distribution  between  zero  and  one. 
The  second  method  systematically  selects  the  ray  direction  to  pass  through 
the  centers  of  equal  solid  angles  which  comprise  the  unit  sphere  or  a 
portion  thereof. 

After  all  ray  directions  are  generated,  tracking  through  the  geometry 
is  commenced  for  the  dose  point  defined.  The  geometry  input  on  SYSDA 
unit  has  been  transformed  to  the  dose  point  as  described  in  the  previous 
section. 

For  each  ray  direction,  the  individual  accumulation  of  thickness  of 
material  penetrated  is  generated.  Prior  to  the  accumulation,  each  path 
length  is  converted  to  a  standard-material  equivalent  path  length.  The 
accumulation  is  made  for  each  direction.  Some  unnecessary  tracking  is 
eliminated  by  comparing  the  octant  of  the  shield  with  the  octant  of  the  ray. 
Other  unnecessary  tracking  is  eliminated  by  hit-or-miss  criteria. 

When  tracking  has  been  completed,  the  normalized  geometry  distribu¬ 
tion  is  computed  and  curve  fit  by  the  least-squares  method.  Doses  can  be 
computed  with  dose  vs.  thickness  curves  for  arbitrary  solar-flare  spectra 
by  integration  of  the  product  of  dose  and  the  derivative  of  solid  angle  as  a 
function  of  thickness  (areal  density).  Also,  the  emergent  fluxes  may  be 
computed  from  incident  energy  spectra  by  integrating  the  product  of  the 
emergent  integral  spectrum  and  the  derivative  of  the  thickness  solid-angle 
distribution  function. 

ESDOSE  is  the  major  subroutine  which  controls  the  elemental  volume 
path-length  computations. 

In  ESDOSE,  the  following  input  data  are  read  for  each  run  of  this 
routine: 

ORD  (I,  J)  for  J  =  1,4  and  I  =  1,6,  where  J  =  1,4  refers  to 

the  point  within  the  Ith  plane  of  the  hexahedron 

In  Figure  13,  the  hexahedron  is  defined  by  six  planes  and  eight  points.  The 
points  are  ordered  counterclockwise  or  in  the  direction  of  a  right-hand  screw, 
directed  outward  from  the  hex  volume.  This  direction  is  the  same  as  the 
surface  normal. 
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The  following  sequence  of  points  is  recommended  for  input  of  the  point 

order: 


ORD  (I,J) 

J 

I 

1 

2 

3 

4 

1 

1 

2 

4 

3 

2 

3 

4 

5 

6 

3 

7 

8 

6 

5 

4 

2 

1 

8 

7 

5 

1 

3 

6 

8 

6 

2 

7 

5 

4 

Next,  the  parameter  IRS  is  read  as  a  dummy  variable.  For  each  areal> 
density  computation,  the  following  input  data  are  read: 

RANDOM  =  indicator  for  ray  selection  method  (>0  for  random, 

=  0  for  equal  solid  angles,  <  0  return  to  the  main  program) 

N  =  number  of  rays  to  be  selected 

NMAT  =  number  of  different  materials 

NRHO  =  number  of  solid  elemental  volumes  whose  density  will  be 
changed 

NRPT1  =  printout  indicator  for  ray  direction  cosines  (x  0,  direction 
cosines  are  to  be  printed;  =  0,  direction  cosines  will  not  be  printed) 

NPRT2  =  printout  indicator  for  tracking  information  (*0,  shield 
identification,  tracking  direction,  equivalent  thickness  and  accumulative 
thickness  will  be  printed;  =  0,  the  previous  quantities  will  not  be 
printed) 

NDIV  =  number  of  microscopic  solid-angle  increments 
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IPUNCH  -  indicator  for  medium  of  sorted  output  (  =  0  printed  only; 

=  7  printed  and  punched;  =  8  printed  and  on  binary  tape;  >0  and  <  7, 
printed,  punched,  and  on  binary  tape) 

If  NRHO  >  0,  the  information  for  altering  the  original  material 
densities  is  read  as  follows: 

[NSS(I)  for  1=1,  NRHO]  where  NSS  is  the  shield  serial  numbers 
lor  which  the  densities  are  changed,  and  [DENS(I)  for  1=1,  NRNO] 
where  DENS  is  new  densities  (gm/cm^)  corresponding  to  the  NSS 
serial  numbers 

CASE  =  alphanumeric  description  of  the  particular  problem  being 
solved 

NEWF9  =  number  of  end  of  file  marks  preceding  the  file  for  writing 
the  areal-density  distribution  function  and  associated  variables  on 
tape  9  (reserve  tape)  {if  NEWF9  =  0,  the  distribution  function  is  not 
saved.  In  order  to  create  the  tape  on  unit  9,  NEWF9  should  be 
input  as  -O 

TMX  =  maximum  areal  density  (gm/cm4)  for  which  tracking  is 
allowed 

EPSLN  =  parameter  used  to  compensate  for  tracking  and  hit-or-miss 
test  round-off  errors 

Next,  the  range-energy  constants  for  the  standard  material  and  the  NMAT 
materials  are  input. 

DSTD  =  8 

s 

ESTD  =  n 

s 

[  DELTA(I),  ETA(I),  I  =  1,  NMAT]  =  [  6A  and  lA  for  NMAT 
materials] 

Titles  for  the  standard  and  actual  materials  are  read  as  [TTLSTD(I), 

I  =  1,3]  and  ( TTLMAT(I,  J),  I  =  1,3  and  J  =  1,  NMAT]. 

The  variable  RANDOM  is  tested  to  determine  whether  random  or 
systematic  selection  of  the  rays  will  be  used. 
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a.  Systematic  Selection 


For  the  systematic  selection  process,  the  ray  direction  cosines  are 
selected  at  the  center  of  equal  solid  angles  with  Equations  4  through  9.  The 
input  angular  boundaries  for  the  azimuthal  (0)  and  colatitude  (6)  angles  are 


input  and  defined  as  fallows : 

PHII 

= 

{angles 

in  deg 

PH  IF 

= 

*F 

THETI 

= 

6I 

THETF 

= 

0F 

The  computed  variables  are 

as 

follows : 

NTSA  = 

NTS 

NTH  = 

Ne 

NPH  = 

N 

0 

THK  (K)  = 

6K 

PHL  (L)  = 

The  systematically  selected  direction  cosines  are  computed  and  defined  as 
follows : 


DIRCOS  (1,J)  = 
DIRCOS  (2,  J)  =  ^ 
DIRCOS  (3,  J)  =  Yj 


where  J  refers  to  the  Jth  ray. 

The  direction  cosines  are  computed  for  NDIV  macroscopic  solid-angle 
increments.  The  NRAYS  direction  cosines  are  printed{when  NPRT1  =  1) 
before  returning  to  ESDOSE. 
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b.  Random  Selection  of  N  Ray  Direction  Cosines 


Subroutine  RNDSEL  is  called  if  RANDOM  >  0.  The  statement 
DUMY  =  RANF(I)  initializes  the  random  number  generation  process. 

The  statement  x(J)  =  RANF(O)  generates  the  random  numbers  where 
0  <  X(J)  <  1  and  J  refers  to  the  three  component  direction  cosines  of  a 

randomly  oriented  vector.  From  Equation  1,  the  unnormalized  direction 
cosines  (oq*  0o»  Vq)  are  computed  as  follows: 

UNDCS(J)  =1.0  -  2.0  *  x(J) 


for  J  =  1,  3. 

The  normalized  ray  direction  cosines  are  defined  as  DIRCOS(J,I)  where 
J  and  I  refer  to  the  Jth  component  of  the  Ith  direction.  After  the  direction 
cosines  are  computed,  they  are  printed  if  NPRT1  =  1.  The  last  computa¬ 
tion  is  DUMY  =  RANF(l),  which  is  used  to  generate  a  new  sequence  of 
random  numbers. 

After  R21,  ESDOSE  calls  the  subroutine  TRACK  to  ray  trace  the 
transformed  geometry  with  the  selected  ray  directions. 

8.  SUBROUTINE  TRACK 

The  octants  of  all  rays  are  computed  by  the  subroutine  OCTCOS  and 
stored  in  the  array  OCTR  for  NRAYS  rays. 

The  variable  IS  is  used  to  count  up  to  the  total  number  of  elemental  vol¬ 
umes  (NOS)  to  be  tracked.  Then  the  first  elemental  volume  parameters  of  a 
shield  are  read  and  stored  in  the  array  SHIELD  (I,  1),  where  the  number  of 
components  of  the  composite  shield  is  defined  as  NCOMP.  Afterwards,  the 
data  for  the  remaining  NCOMP  -  1  volumes  are  read  and  stored  in  the  array 
SHIELD  (I, ICP).  The  shield  serial  numbers  are  compared  with  NSS(I)  if 
NRHO  >  0  and  densities  are  changed  as  required. 

The  tracking  process  begins  by  selecting  ray  number  IR.  The  array 
PATH  for  NRAYS  will  contain  the  accumulative  areal  density.  If  PATH(IR) 

2:  the  prespecified  maximum  areal  density,  tracking  is  discontinued  for  the 
IRth  ray. 

Then  variables  IPLS,  ING,  and  INTS  are  initialized  as  zero,  and  used 
respectively  to  count  the  number  of  positive,  negative,  and  total  number  of 
elemental  volumes  intersected  for  the  IRth  ray. 
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The  ray  direction  cosines  are  stored  in  SHLDD  (38-40)  for  all  volume 
path-length  computation  subroutines. 

Then,  R22  tracks  through  each  composite  shield  component  for 
the  IRth  ray.  The  elemental  volume' data  SHIELD  (I,  ICP)  is  then  transferred 
to  SHLDD  (I).  Density  (RHO),  material  type  =  [  SHLDD  (2)  +  l],and 

shield  octant  (SOCT)  are  then  defined.  If  SOCT  =  9,  the  ray  tracing 
procedure  will  continue.  Otherwise  SOCT  is  compared  to  the  ray  octant 
[OCTR  (IR)  ] .  If  they  are  not  equivelent,  the  tracking  procedure  is  discon¬ 
tinued  and  another  volume  is  selected  for  ray  tracing. 

When  the  tracking  process  continues,  the  volume  path  length  (SHLP) 
is  set  to  zero  before  the  volume  path-length  computation  begins.  Then  CS  is 
defined  as  ±  1.0  to  designate  a  solid  and  a  void  region,  respectively.  After 
R23,  a  logical  test  is  used  to  eliminate  unnecessary  tracking.  If  IPLS  =  0 
and  CS  =  -1,  there  are  no  positive  composite  shield  components.  The 
volume  path  length  is  now  computed  by  one  of  the  following  subroutines: 


Volume 

Subroutine 

Hexahedron 

TKHEX 

Cylinder 

TKCYL 

Sphere 

TKSPH 

Hemisphere 

TKHM 

Cone 

TKCON 

Truncated  cone 

TKCON 

Ellipsoid 

TKELL 

The  arguments  of  the  subroutines  are  defined  as  follows: 

SHLP  =  path  length  computed  through  the  volume 

DP(1),  DP(2)  =  minimum  and  maximum  distances  through  the  volume 
for  the  IRth  ray 

Hit-or-miss  logical  tests  are  also  made  in  each  subroutine  to  avoid 
unnecessary  tracking.  If  SHLP  is  zero  after  R24,  etc.,  control  is  trans¬ 
ferred  to  R25  because  there  is  no  volume  intersection,  and  the  next  volume 
is  selected.  If  not  zero,  the  distances  through  the  volume  surfaces  are 
converted  from  inches  to  centimeters,  designated  as  plus  or  minus,  and 
stored  in  the  array  ELEMP.  The  plus  and  minus  signs  refer  to  solid  and 
void  shields  respectively.  RHO  and  MAT  are  stored  as  RHOS  and  MATS. 
Then  INTS  and  ING  or  IPLS  are  incremented.  Next,  the  subroutine  ELIM? 
is  called  to  determine  if  any  or  all  negative  volumes  completely  void  all  the 
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solid  regions.  If  the  solid  regions  of  a  composite  are  voided,  the  next 
ray  is  selected  and  ray  tracing  is  continued.  Otherwise,  the  areal  density 
(XS)  is  computed  for  the  composite  shield  components  by  subroutine 
COMPSP.  The  areal  density  is  then  accumulated  in  the  array  designated 
as  PATH. 

The  above  process  continues  until  computations  for  all  NRAYS  have 
been  completed. 

Then  IS,  which  is  now  the  total  number  of  elemental  volumes  processed, 
is  increased  by  (NCOMP-1).  If  IS  is  less  than  NSO  (the  total  number  of 
volumes),  another  composite  shield  volume  set  is  selected  for  tracking. 
Otherwise  control  is  returned  to  subroutine  ESDOSE. 

The  subroutines  used  by  TRACK  will  now  be  explained  in  more 

detail. 

9.  SUBROUTINE  OCTCOS 

This  subroutine  uses  the  same  method  described  in  Section  III-4.b.  to 
compute  octant  numbers.  However,  the  previous  computations  were  for 
volume  rather  than  for  ray  octants.  For  zero  radii  the  two  methods  are 
the  same. 

10.  SUBROUTINE  ELIMS 

After  the  solid  shields  of  a  composite  shield  have  been  tracked,  the 
routine  ELIMS  is  used  after  tracking  each  void  component  to  determine  if 
the  solid  shield  traversal  paths  have  been  voided.  This  eliminates  unneces¬ 
sary  ray  tracing  and  areal-density  computations  for  some  of  the  composite 
shield's  void  components.  The  array  ELEMP(I,  ICP)  contains  the  entrance 
(1=1)  and  exit  (1=2)  distances  to  ICP  composite  shield  components.  The 
array  elements  are  positive  and  negative  for  solid  and  void  regions  respec¬ 
tively.  Logic  used  in  ELIMS  is  discussed  in  detail  in  Section  II- 6 .  c.  The 
variable  ISIT  is  initially  zero  and  changed  to  unity  if  the  solid  shield  com¬ 
ponents  are  voided  (R26). 

In  the  loop  R27  the  distance  for  the  solid  (positive)  shields  are 
transferred  from  ELEMP  to  DPPOS.  ING,  the  number  of  void  components 
that  have  been  tracked,  is  tested  for  a  nonzero  value.  If  solid  and  void 
components  have  been  tracked,  R28  orders  the  void  region  distances  by 
entrance  path  lengths  in  array  DPNEG.  Then  R29  independently  deter¬ 
mines  the  remaining  portions  of  each  solid  shield  component. 
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11.  SUBROUTINE  COMPSP 


If  the  solid  portions  of  a  composite  shield  arc  not  voided,  subroutine 
TRACK  (R30)  calls  COMPSP  to  compute  the  solid  path  areal  density. 

COMPSP  is  also  called  for  single-element  solid  shields.  This  routine 
computes  the  following:  the  distance  from  the  dosimeter  to  the  exit  surface 
of  each  solid  region  (X(l)  or  X(3)];  the  accumulative  standard-material  areal 
density  [XS],  and  the  solid  region's  areal  densities  [DM]  for  each  solid  shield 
component.  These  data  are  stored  on  SYSDA  20. 

12.  SUBROUTINE  TKHEX  (HXPTH,  NP,  DP) 

NP  and  DP  have  been  defined  as  the  order  of  the  hexahedron  coordinate 
points  and  the  path  lengths  to  the  hexahedron  surfaces.  HXPTH,  the  path 
length  through  the  hexahedron  volume,  is 

HXPTH  =  jDP(2)  -  DP(1)| 

If  the  dosimeter  is  inside  the  hexahedron,  DP(1)  is  defined  as  zero  and  DP(2) 
is  the  distance  from  the  dosimeter  to  surface. 

NT  specifies  whether  the  subroutine  logic  is  used  to  determine  if  the 
dosimeter  is  inside  (NT=2)  or  to  compute  DP  (NT=1). 

NT  is  set  equal  to  unity,  and  the  number  of  non-zero  positive  intersec¬ 
tions  with  the  hexahedron  faces  is  calculated.  If  the  number  of  intersections 
is  unity,  the  dosimeter  is  inside  and  HXPTH=DP(2).  The  hit-or-miss 
criteria  eliminate  cases  with  no  intersections.  If  there  are  two  or  more 
intersections,  there  are  possible  path-length  equivalent  combinations  with 
different  volume  path  lengths.  For  example,  there  may  be  glancing  incidence 
at  a  point  on  the  surface  or  a  ray  from  a  dosimeter  inside  the  hexahedron 
volume  passing  through  an  apex  or  corner  which  has  the  same  surface  path 
lengths.  These  ambiguities  will  be  eliminated  by  determining  if  the  dosimeter 
is  inside  or  outside  the  hex  volumes,  and  using  the  following  logic: 


Type  of  Intersection 

Volume  Path  Length 

Dosimeter  inside  and  passing  through  a 
corner  or  apex  point 

HXPTH  =  DP(2) 

Dosimeter  outside  with  glancing 
incidence 

HXPTH  =  DP(2)  -  DP(1)  =  0 

Dosimeter  outside  with  ray  in  the  plane 
of  a  hexahedron  face 

HXPTH  =  DP(2)  -  DP(1)  /  0 
(see  the  second  paragraph  of 
subsection  II-6.b. ) 
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The  subroutine  TKHEX  start®  by  initializing  XO  and  defining 
DIR  =  DIRCSA. 

XO(K),  K=l,  3  are  the  coordinate  components  of  point  M  of  Figure  13. 

DIRCSA(K)  are  the  direction  cosines  of  the  tracking  ray. 

NT  is  set  equal  to  unity  to  indicate  that  the  surface  path  lengths  will 
be  computed  for  the  tracking  ray.  Equations  35,  36  are  then  used  to  com¬ 
pute  the  trigonometric  functions  for  the  rotation  matrix  M,  (Equation  34).  In 
the  program 


A(I,J)  -  M(I,J) 

At  R31  the  DSCS  coordinates  X  (I,  J)  are  transformed  to  the  RDCS  as 
XR(I,  J).  I  and  J  refer  to  the  coordinate  point  and  component  respectively. 
Next  NT  is  tested,  because  if  the  calculation  mode  is  determining  whether 
the  dosimeter  is  inside  or  outside  the  hexahedron  volume,  hit-or-miss 
criteria  are  not  necessary.  If  NT  =  2  it  has  been  determined  that  the  ray 
intersects  the  volume.  If  it  had  not  intersected  the  volume  there  would  have 
been  an  exit  near  R32  ,  Between  R33  and  R34  the  hexahedron  miss -rectangle 
criteria  (Figure  16)  are  used  in  the  form  of  the  following  logic  and  definitions: 

Let  x(J,  X)  be  defined  as  the  Kth  component  of  the  Jth  point  of  a  pro¬ 
jected  hexahedron  face  in  the  RDCS  (x,y)  plane 

If  x(J,  K)  is  bounded  by  ±  *  ,  the  coordinate  x  of  point  J  is  at  a 
neutral  position  and  is  considered  both  negative  and  positive 

The  variables  NEG  and  IMPOS  are  incremented  for  negative  and  posi¬ 
tive  coordinate  components,  respectively.  If  all  XR(J, K),  for  a  given  K, 
have  the  same  value,  the  miss  rectangle  will  not  contain  the  origin,  i.  e.  , 
the  ray  misses  the  rectangle.  The  two  tests  between  R33  and  R34  check 
for  this  condition.  If 


IMPOS  or  NEG  =  8 

for  either  the  XR  or  YR  component,  the  ray  misses  the  rectangle  and  control 
is  transferred  back  to  TRACK.  Otherwise,  the  calculation  is  continued. 

INTSCT  as  the  number  of  intersections  with  hexahedron  planes 

NPL  *  number  of  one  of  6  hexahedron  faces 

Now  the  hexahedron  miss  rectangle  will  be  applied  to  each  face  of  the  hexa¬ 
hedron.  If  there  is  a  miss  for  the  NPL  plane,  i.  e.  ,  IMPOS  or  NEG  =  4, 
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there  is  a  transfer  to  R35  and  the  next  face  is  tested  for  intersection.  If 
there  is  a  possible  hit,  i.  e.  ,  IMPOS  and  NEG  <  4,  transfer  is  made  to  R36  , 
which  makes  the  final  miss  test.  This  test  (Equation  51)  performs  an  exact 
test  for  determining  whether  there  is  a  hit  or  miss.  The  following  definitions 
are  used  to  implement  this  test: 

Jl,  J2  =  points  defining  the  vectors  VR(J1)  and  VR(J2) 

VECTPR  =  IT  •  [vR (Jl)  x  VR(J2)j  (Figure  11) 

After  R37,  IMPOS  and  NEG  are  tested  for  an  intersections  and  INTSCT  is 
incremented  if  the  test  is  satisfied.  Also  the  array  NFACE  is  defined  such 
that  NFACE(I)  is  equal  to  the  set  of  planes  intersected,  with  1=1,  INTSCT. 

If  NT=1  ,  this  is  the  first  cycle  through  the  subroutine  and  INTSCT  is 
checked  for  plane  intersections.  If  INTSCT=0  ,  there  are  no  intersections 
and  there  is  a  return  to  TRACK.  Otherwise,  R38  computes  the  path  lengths 
to  the  intersected  planes.  R39  redefines  the  hexahedron  points  in  the  DSCS 
as  XC(J,  M).  Then  R40  computes  the  components  of  vectors  (Figure  10) 


A  „  nd  B,  where 


A  =  V(l)  -  V(2) 
B  =  V(3)  -  V(2) 


and 


the  Jth  component  of  V(I)  =  XC(I,  J) 

The  vector  AxB  is  in  the  direction  of  the  plane  defined  by  AxB,  so  ep^  of 
Equation  23  is  computed  from 

-  AxB 

e _ =  : - 

where 

AKB(I)  = 

VAKB  S 


PL, 


I  AxB  I 


Ith  component  of  AxB 
IaxbI 
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The  components  of  epL  (R41),  i.e.,  the  direction  cosines  of  the  hexahedron 
planes, are  defined  as  ALPHA(M),so 

a  =  ALPHA(l),  etc. 


From  Equation  23 


a  a  +SB  +  v  V 
L  PL  HL  PPL  YL  PL 


DEN 


and 


ALPHA(I)  =  Ith  component  of  e*! 

If  |DEN|  <  c  where  «  is  =  10“^,  i.  e.  ,  DEN  =  0,  the  line  e".  is  parallel  to 
the  hexahedron  plane  and  there  is  no  intersection.  In  this  case  the  path 
length  (PtH  =  PR)  is  set  equal  to  a  large  number,  namely  10^®.  Later,  such 
large  path  lengths  will  indicate  plane  misses  or  non-intersections.  This 
large  value  will  also  be  substituted  for  negative  path  lengths  or  PR(I)  <  0. 

The  path  length  (PR)  was  computed  (R43)  from  the  following  form  of 
Equation  49: 


P 


TH 


since 


=  0 


After  R42  if  NT  =  1  ,  INTSCT  is  checked  for  more  than  one  inter¬ 
section.  If  there  is  only  one  intersection  it  is  equivalent  to  the  volume  path 
length  because  the  dosimeter  must  be  inside  the  hexahedron.  For  several 
intersections,  R44  uses  PMIN  =  10^0  and  the  library  routine  AMIN1  is 
used  with  PR(I)  to  determine  the  minimum  path  length  PMIN  =  DP(1).  It 
also  stores  the  array  PR(I)  into  PRHEX(I). 

R45  replaces  the  minimum  PR(I)  by  10^0.  R46  computes  the 

coordinates  of  the  point  M  (Figure  13).  DIST  is  the  distance  from  the  dosim¬ 
eter  to  M.  The  R47  statements  compute  the  direction  cosines  of  a  vector 
from  M  to  D  and  translate  the  origin  to  point  M.  NT  has  been  set  equal  to  two 
and  control  is  now  transferred  to  R48.  This  recycle  will  compute  the  allowed 
path  lengths  to  the  hexahedron  surfaces  for  a  ray  directed  from  M  to  D. 

The  minimum  path  length  can  be  compared  with  DIST  (the  distance  from  M 
to  D)  to  determine  whether  the  dosimeter  is  inside  or  outside  the  hexahedron 
volume.  The  cal«.  Nations  will  proceed  as  before  up  to  R43.  Then  NT  will 
be  tested  and  con..ol  transferred  to  R49  ,  which  will  compute  the  minimum 


path  for  the  vector  directed  from  M  to  D.  If  MP  >  MD,  the  dosimeter  is 
inside  and  the  volume  path  length  HXPTH  is  PMIN.  In  this  case  the  ray 
passes  through  a  corner  or  edge  of  the  hexahedron.  Otherwise,  K50 
computes  the  second  minimum  and  the  volume  path  length  is  the  difference 
between  the  first  and  second  minima.  If  there  is  glancing  incident  e,  the 
two  path-length  minima  will  be  identical  and  the  volume  path  equals  zero. 

13.  SUBROUTINE  TKCYL 

Since  tracking  through  all  shields,  except  hexahedrons  and  spheres, 
is  done  in  the  RDCS,  the  ray  direction  cosines  DIRSCA  are  rotated  by  the 
rotation  matrix  A(K,  J)  for  cylinder s, cones,  and  ellipsoids,  etc.  ,  as 
follows: 


DIRRAY(K) 


s 


L=1 


A(K,  L)  DIRCSA(L) 


In  the  rotated  system,  the  z  axis  of  the  cylinder  is  parallel  to  the  coordinate 
z  axis.  Starting  at  R51  the  calculations  will  determine  whether  the  dosim¬ 
eter  is  inside  or  outside  the  cylinder.  If  it  is  inside  (Figure  8),  the  origin 
(dosimeter  locatior  )  is  bounded  by  the  z  component  of  points  1  and  2.  In 
the  RDCS  the  cylincer  coordinates  are  defined  by  XCY(I,  J),where  I  and  J 
refer  respectively  to  the  point  and  component.  If  the  dosimeter  is  inside, 

C  of  Equation  21  is  negative  or  zero  (R52).  If  the  detector  is  inside  the 
cylinder,  the  path-length  computation  is  as  follows: 

If  DIRRAY(3)  =  0,  the  ray  is  in  the  (x  ,  y)  plane  and  therefore  will 
not  intersect  either  of  the  plane  bases.  In  this  case  control  is  transferred 
to  R5^  and  Equation  47  is  implemented.  The  components  for  the  quadratic 
solution  for  intersection  with  the  cylinder  are  evaluated  with  subroutine 
ROOT.  The  subroutine  ROOT  solves  the  conic  surface  quadratic  for  the 
path  lengths  =  PR(3),  PR(4).  If  a  path  length  is  negative,  it  is  defined  as 
1020,  since  it  is  not  allowed.  The  minimum  of  the  roots  PR(3),  PR(4)  is 
the  path  length.  If  DIRRAY(3)  *  0,  the  ray  intersects  the  cylindrical  surface 
or  a  base  plane.  The  computation  is  as  follows:  for  each  base  1=1,  2  the 
distance  from  the  detector  to  the  bases  are  computed  from  Equation  40, 
whose  computer  equivalent  is 


DP(I) 


ZBOUND(I) 

DIRRAY(3) 
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if  DP(I)  <  0,  the  ray  does  not  intersect  the  Ith  base.  Otherwise,  since  the 
dosimeter  is  inside,  the  projection  of  the  ray  on  the  (x  ,  y)  plane  must  be 
less  than  or  equal  to  the  cylinder  radius  for  intersection.  The  end  point 
coordinates  (R54)  of  the  projected  line  are  defined  as  CXY(J),  and  C  of 
Equation  2.1  is  defined  as  C2.  If  C2  is  positive,  the  ray  intersects  the 
cylindrical  surface  and  control  is  transferred  to  R53.  The  subsequent 
operations  have  been  explained.  If  the  detector  is  outside,  all  possible 
types  of  intersection  must  be  examined.  If  DIRRAY(3)  =  0,  intersection 
possibilities  exist  only  with  the  cylindrical  surface.  Therefore,  path  lengths 
PR(I),  1=1,  2  are  defined  as  10^0.  if  DIRRAY(3)  *  0,  intersection  possi¬ 
bilities  with  the  cylinder  bases  are  considered.  The  distances  to  the  bases 
are  computed  as  before.  From  Figure  21  the  projected  center  and  point  of 
intersection  of  the  ray  with  the  base  must  be  separated  by  a  distance  less 
than  the  cylinder  radius,  or  the  associated  path  length  is  not  allowed,  i.  e.  , 
Rdc»  which  is  defined  as  C2,  is  less  than  R.  The  next  section  (R55) 
computes  the  path  lengths  for  intersection  with  the  cylindrical  surface.  If 
YR  of  the  tracking  ray  is  1.  0  or  A  of  Equation  21  is  zero,  there  is  no  inter¬ 
section.  Otherwise  the  quadratic  solution  must  be  used  to  determine  the 
cylinder  path  lengths  PR(3)  and  PR(4).  The  coordinates  of  the  points  of 
intersection  of  the  tracking  ray  with  the  cylinder  (R56)  are  computed  and 
checked  to  determine  if  they  are  bounded  by  the  coordinates  of  the  cylinder 
end  surfaces.  If  the  z  coordinates  are  not  bounded,  the  corresponding  path 
lengths  are  not  allowed.  The  two  minima  of  the  PR(I),  1=1,  4  are  computed 
(R57)  and  the  path  (CYPATH)  through  the  cylinder  is  the  positive  difference 
between  them. 

14.  SUBROUTINE  TKSPH 

The  quadratic  form  of  solution  is  evaluated  in  the  DSCS  and  invalid 
roots  are  set  to  the  arbitrarily  large  value.  The  distance  to  the  center  of 
the  sphere  is  compared  to  the  sphere  radius  to  determine  the  detector  loca¬ 
tion.  If  inside,  the  path  length  is  equivalent  to  the  minimum  value  computed 
by  the  root  evaluation  (Equation  18).  If  outside,  the  volume  path  length  is  the 
positive  difference  between  the  two  spherical  surface  path  lengths. 

15.  SUBROUTINE  TKHM 

The  ray  direction  cosines  are  transformed  to  the  RDCS(R58)  and  path 
lengths  to  the  spherical  surface  are  computed  by  ROOT.  ISC  =  0,  1, 
respectively,indicates  that  the  ray  hits  or  misses  the  spherical  surface. 

If  ISC  =  1,  the  ray  must  miss  the  hemisphere  which  is  contained  by  the 
sphere.  If  the  spherical  surface  paths  are  valid  (ISC  =  0),  R59  computes  the 
coordinates  of  the  points  of  intersection  with  the  spherical  surface.  Between 
R59  and  R60  tests  are  made  to  determine  if  the  intersection  coordinates  are 
bounded  by  points  1  and  2  (Figure  8).  Next  the  hemisphere  base  is  examined 
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Figure  21.  Projection  of  Cylinder  Base  Plane 


for  possible  intersection  as  follows:  if  DIRRAY(3)  =  0,  there  is  no  inter¬ 
section  with  the  base  and  PR(2)  is  set  equal  to  10^0. 

If  DIRRAY(3)  *  0,.  the  distance  R61  to  the  plane  DPL  is  evaluated. 

If  DPL  is  within  the  radial  bounds  and  positive,  it  is  a  valid  path  (R62),  and 
the  two  minima  (not  equal)  of  PR(2),  PR(3),  and  PR(4)  are  selected  as 
BP(1)  and  DP(2).  From  Figure  8,  if  the  origin  is  bounded  by  the  z  (R63) 
components  of  points  one  and  two,  the  dosimeter  is  inside  and  the  hemi¬ 
sphere  volume  path  length  is  DP(1).  If  the  detector  is  outside  and  DP(1)  is 
10^0,  there  is  no  intersection.  Otherwise  the  path  length  is  the  difference 
between  the  two  minima. 

16.  SUBROUTINE  TKCON 

NP  is  used  to  designate  the  type  of  cone.  Values  of  2  and  3  refer 
respectively  to  non-truncated  and  truncated  cones.  R64  defines  the 
cone's  RDCS  coordinates  X  (I,  J).  R65  defines  the  rotation  matrix  AR 

(I,J)  which  is  used  in  R66  to  transform  the  ray  direction  cosines  from  the 
DSCS  to  the  RDCS.  Also  the  cone  radii  at  the  base  and  truncation  plane  are 
respectively  defined  as  RB  and  RT.  After  R67,  logical  tests  are  used  to 
determine  whether  the  dosimeter  is  inside  or  outside  the  cone  volume.  If 
the  dosimeter  is  inside  (Equation  40)  the  volume,  the  origin  is  bounded  by 
the  z  components  of  Zr(1)  and  Zr(NP)  and  the  distance  from  the  center  line 
of  the  cone  to  the  dosimeter  or  origin  must  be  less  than  the  radius  of  the 
cone  cross  section  in  the  (xj^,  y^)  plane.; 

RCD  =  radius  of  the  cone  cross  section  in  the  (x,  y)  plane 


/RDET  =  distance  from  the  cone  centerline  to  the  dosimeter 
If  the  detector  is  inside,  special  cases  of  Yr  are  considered  as  follows: 

Yj^  =  0,  the  ray  is  in  the  (xr,  y^)  plane  and  only  intersects  the 
conical  surface  (R68) 

Y_  =  -1,  the  ray  intersects  the  base  plane  (R69) 

R 

and  for  the  more  general  case  the  path  lengths  to  the  bases  (Zl)  are  computed 
first  (R70).  Then  they  are  tested  for  positive  values. 

At  R71  the  conical  surface  path  lengths  [PR(3),  PR(4)]  are  computed 
and  the  minimum  (R72)  of  the  three  path  lengths  is  selected  as  the  volume 
path  length.  If  the  dosimeter  is  outside  the  cone  volume,  control  is  trans¬ 
ferred  to  R73.  Then  Yr  is  tested  for  a  zero  value.  If  zero,  the  ray  is  in 
the  (xr,  ypj  plane  and  does  not  intersect  either  of  the  cone  bases.  Other¬ 
wise,  the  distance  to  the  base  plane  PR(1)  is  computed  and  tested  for  a 
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negative  value.  Then  PROJR  must  be  less  than  RB^.  PROJR  is  the  square 
of  the  (xr,  y^)  plane-projected  distance  of  a  vector  from  the  plane  point  of 
intersection  to  the  cone  centerline.  A  similar  method  is  used  to  compute 
the  path  length  PR(2)  to  the  truncation  plane.  Subsequently  at  R74  the  path 
lengths  to  the  conical  surface  PR{3),  PR(4)  are  calculated.  Starting  at 
R75  ,  the  two  minimum  path  lengths  are  computed  and  the  volume  path 
length  defined  as  their  difference. 

17.  SUBROUTINE  TKELL 

At  R76  ,  the  ray  direction  cosines  of  the  ellipsoid  coordinate  axes 
are  transformed  by  DIR  (I,  J)  to  the  RDCS  as  DIRRAY(I).  The  transforma¬ 
tion  matrix  is  defined  by  Equation  31.  JMAX  is  used  to  indicate  the  type  of 
ellipsoid  and  after  R77  the  z  bounds  of  the  ellipsoid  are  established  in  the 
RDCS.  JMAX  and  the  bounds  are  defined  as  follows  (Figure  6): 

JMAX  =  4,  the  ellipsoid  is  nontruncated 

=  5,  the  ellipsoid  is  truncatt  i  and  bounded  by  the  coordinates  of 
points  3  and  5 

=  -5,  truncated  and  bounded  by  point  5  and  ZB(1)  =  X(4)  -  C 

=  6,  truncated  and  bounded  by  points  5  and  6 

After  the  z  coordinate  bounds  are  specified.  R78  computes  the  path  lengths 
to  intersections  with  the  coordinate  planes. 

If  IEL  =1,  2,  3,  which  is  equivalent  to  |jMAX|  =  4,  5,  6,  (R79)  there 
are  0,  1,  and  2  intersections  with  the  truncation  planes.  If  there  are  two 
plane  intersections,  and[PR(l)  +  PR(2)J  >  10^0  there  is  at  most  one  inter¬ 
section  with  the  truncation  planes.  In  this  case,  a  transfer  is  made  to  R80 
to  find  a  possible  intersection  with  the  ellipsoid  surface.  This  subroutine 
allows  for  cases  where  the  dosimeter  is  located  inside  the  ellipsoid. 

At  R81  ,  the  z  coordinate  of  the  ellipsoid  intersections  are  checked 
to  see  if  they  are  bounded  by  the  z  coordinates  of  the  truncation  planes. 

R82  selects  two  nonequal  minimum  path  lengths,  i.  e. ,  if  such  exist. 

At  R83,  C  of  Equation  17  is  checked  for  a  negative  value.  If  C  is  negative, 
the  dosimeter  is  inside  the  ellipse  and  the  volume  path  length  is  the  first 
minimum  of  the  path  lengths.  If  the  dosimeter  is  outside  the  ellipsoid  and 
either  of  the  two  minimum  path  lengths  is  10^®,  the  ray  misses  the  ellipsoid. 
Otherwise,  the  volume  path  length  is  the  difference  between  two  minima. 

At  R84  of  subroutine  TRACK  a  transfer  is  made  to  subroutine  GEOMDS 
to  compute  the  normalized  geometry  distribution  and  curve  fit  coefficients. 
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18.  SUBROUTINE  GEOMDS 


In  order  to  represent  the  areal-density  distribution  function  more 
accuratelyt  least-square  curve  fits  are  made  by  fitting  the  function  in  sections 
for  particular  thickness  intervals.  This  is  necessary  because  a  large  por¬ 
tion  of  the  primary  dose  comes  from  the  region  of  smaller  thicknesses  and 
normally  there  is  a  large  variation  in  the  distribution  from  minimum  to 
maximum  thickness. 

The  input  data  variables  are  defined  as  follows: 

NSECT  =  number  of  curve  fit  sections 
NSDEG(N)  =  the  degree  of  polynomial  used  for  each  section 
XMAX  =  TMX  (gm/cm2) 

TAU  =  t 

DELR  =  6R 

TMAX  =  XMAX  =  maximum  areal  density  for  least-squares  curve 

fits 


The  minimum  value  of  THK,  TMIN  is  selected  from  the  computed  equivalent 
thicknesses  previously  generated.  TMIN  is  used  with  TMAX,  t,  and  6R  to 
compute  the  weighted  THK  values  with  Equation  61. 


The  maximum  thicknesses  (XFIT)  are  required  for  each  section  to  be 
curve  fit  and  are  computed  using  the  logarithmic  relationship 


XFIT(L) 


TMIN 


/TMAX  \ 
\  TMIN  / 


L/ NSECT 


R85  computes  the  fraction  of  the  solid  angles  (GM)  with  thickness  less 
than  THK. 

The  normalized  distribution  is  computed  by  dividing  the  GM  array  by 
the  number  of  rays  (NSA)  in  R86  .  .. 

The  loop  R87  computes  the  least-square  polynomial  curve  fit 
coefficients  (CSECT)  for  each  areal-density  section.  A  classical  least- 
squares  subroutine!  designated  LEAST,  is  used. 
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19.  SUBROUTINE  ORDER 


Subroutine  ORDER  sorts  the  path-length  data,  which  have  been  stored 
on  SYSDA  unit  20  by  subroutine  COMPSP,  according  to  ray  number  and 
decreasing  distance  relative  to  dosimeter  location.  The  following  variable 
definitions  are  used: 


NMAX  =  the  maximum  number  of  core  storage  locations 
available  for  the  sorted  data 


BPL  =  the  array  which  contains  the  sorted  distances,  path 
lengths  and  material  numbers  data;  the  data  are 
packed  to  use  all  available  storage  locations  and  are 
dimensioned  by  NMAX  (21000) 

IPATH  =  total  number  of  path  lengths  computed  for  each  ray 
(dimensioned  1000) 

IR1  and  IR2  =  the  minimum  and  maximum  ray  numbers  for  which 

core  is  available  for  storing  path  length  data  into  BPL} 
these  parameters  may  not  necessarily  be  1  and 
NRAYS,  respectively,  since  the  number  of  locations 
needed  for  the  complete  set  of  data  can  exceed  NMAX 


The  first  section  ORDER  performs  initialization  of  IR1  as  1,  zeros 
out  arrays  BPL,  PLDATA  and  MATNO  (final  sorted  data  per  ray)  and 
computes  the  number  of  unsorted  data  records  (NEND)  on  the  temporary 
storage  SYSDA  unit  20.  The  number  of  records  (R88)  is  computed  as 

NRAYS 

NEND  =  J  IPATH(I) 
i=l 


The  next  sequence  of  operations  (R89-R90)  determines  the  number  of  storage 
locations  (NSMAX)  used  in  BPL  for  the  data  for  rays  numbered  IR1  to  IR2. 
The  NSMAX  summation  is  continued  until  NSMAX  exceeds  NMAX  or  IR2 
becomes  equal  to  NRAYS.  Then  the  locations  (IPOSN)  for  storing  the  path- 
length  data  into  the  BPL  array  are  computed. 
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The  loop  R91  reads,  prints,  and  stores  the  path-length  data 
corresponding  to  all  rays  and  paths  contained  in  the  ray  number  interval 
(IR1,  IR2).  These  variables  are  defined  as: 


IR  =  particular  ray  number 
IPT  =  path  number 

X  =  distance  to  the  volume  exit  surface 
DX  =  path  length  (gm/cm^)  through  the  shield 


MAT  =  material  number 


The  next  section  (R92)  performs  the  sorting  according  to  the 
variable  X  for  each  ray  and  produces  the  desired  output  on  the  required 
output  data  medium.  The  sorting  is  executed  as  follows: 

The  number  of  solid  parts  encountered  for  IRth  ray  is  NHIT,  which 
is  picked  up  as  IPATH  (IR).  The  set  of  data  'or  each  ray  is  selected  sequen¬ 
tially  from  BPL.  At  R93  ,  the  sorted  data  results  are  printed  and  stored 
on  the  designated  medium.  The  variable  IPUNCH  determines  the  output  data 


as  follows: 

I  Punch 

Output  Data  Form 

Positive  and  less  than  7 

Punched  cards  on  binary  tape  on  unit  8 

7 

Punched  cards 

8 

Binary  tape  on  unit  8 
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SECTION  IV 


PROGRAM  UTILIZATION  (USER'S  MANUAL) 


This  section  defines  the  procedure  for  data  preparation,  the  MEVDP 
input  data  format,  and  tape  unit  use;  it  is  not  intended  to  be  used  entirely 
independently  of  other  sections  of  this  document. 

1.  NUMERICAL  GEOMETRICAL  DATA  PREPARATION 


A  numerical  geometrical  model  contains  two  types  of  data:  (1)  analyt¬ 
ical  geometrical  coordinates,  and  (2)  material  parameter  data.  The  material 
data  are  the  shield's  material  type  or  code  (MATCD)  and  its  density  (RHO) 
in  units  of  grams /centimeter^.  The  material  code  is  a  number  referring  to 
a  particular  material  type.  The  analytical  coordinate  and  parameter  data 
include  the  type  of  geometrical  volume  (ST),  point  coordinates  in  inches,  and 
volume  radius  (RPHI)  or  cone  half-angle  in  inches  and  radians,  respectively. 

If  the  volume  is  an  ellipsoid,  an  additional  parameter  (JMAX)  defines  its 
degree  of  truncation.  The  elemental  volumes  and  their  defining  points  are 
depicted  in  Figure  4.  AFWL  has  found  it  efficient  to  use  a  standard  form  to 
record  the  data  for  each  shield.  Figure  22  is  a  reproduction  of  the  form 
used  for  this  purpose,  the  Sectoring  Code  Work  Sheet.  Figure  23  is  the 
instructions  printed  on  each  work  sheet  to  assist  in  its  preparation.  The 
data  on  each  work  sheet  are  sufficient  to  define  the  location,  nature,  material, 
and  type  of  shield.  The  relationships  between  the  data  blocks  in  Figure  22 
and  the  computer  code  variables  are  summarized  as  follows: 

SSN  =  shield  serial  number 


ST  = 
MATCD  = 

RPHI  = 
RHO  = 
JMAX  = 


shield  type 

material  code  number  which  defines  the  shield  material 
composition 

radius  of  cylinder  or  radian  half-angle  of  a  cone 
shield  density  (grams/centimeter^) 
type  of  ellipsoid  (degree  of  truncation) 
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SECTORING  CODE  WORK  SHEET 


} 

I 

\ 


SHIELD  SERIAL  NO. 

1 

SHIELD  TYPE 

m 

. 

MATERIAL  CODE 

NO. 

CHEMICAL 

COMPOSITION 

OCTANT  NO. 

DENSITY  (GM/CM3  ) 

CYL  RAD  (in.)  OR  CONE  HALF-ANG  (RAD) 

SPHERE  (2) 

1  X 

Y 

Z 

HEMISPHERE  (2) 

CYLINDER  (2) 

2 

CONE  (2) 

■ 

TRUNCATED 

CONE  (3) 
(TRUNCATION  PT) 

3 

■ 

SIMPLE 

ELLIPSOID  (4) 

4 

1 

TRUNCATED 
ELLIPSOID 
(5  OR  6) 

(TRUNCATION  PT) 

5 

6 

1 

HEXAHEPRON  (8) 

7 

8 

REMARKS 

Figure  22.  Sectoring  Code  Work  Sheet 
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INSTRUCTIONS  FOR  FILLING  OUT  SECTORING  CODE  WORK  SHEETS 


A  right-handed  coordinate  system  is  used,  about  whose  origin  all  shields  are 
described.  All  shields  must  be  one  of  the  geometrical  shapes  listed  below.; 
Coordinates  must  be  in  inches,  angles  in  radians,  and  densities  in  g/crn  . 

1.  Shield  Serial  Number:  use  a  unique  4 -digit  number  for  each  shield.  Numbering 
shall  reflect  a  logical  breakdown  of  the  geometry  being  sectored, i.e.  digits  1 
and  2  denote  the  major  segment,  digit  3  the  order  of  shields  within  the  segment, 
and  digit  4  the  number  of  a  shield  within  a  composite  shield. 

2.  Shield  Type:  ?  2-digit  integer  (see  list  belcw)  denoting  the  shield  type. 

3.  Octant  Number:  a  1-digit  integer  denoting  the  octant  in  which  the  shield  ap¬ 
pears;  if  a  shield  overlaps  more  than  1  octant,  use  octant  number  9. 

4.  Material  Code;  a  2 -digit  number  assigned  to  denote  the  material  type.  The  chem¬ 
ical  composition  shall  be  placed  in  the  space  provided,  if  space  permits,  other¬ 
wise  use  the  remarks  block  (at  bottom  of  form) . 


5-  Cvc Under  Radius. Cone  Half-Angle, Density; 

6.  Shield  Input  Requirements: 

a.  Sphere  (2  points) 

1.  center  of  sphere 

2.  any  point  on  surface 

b.  Hemisphere  (2  points) 

1.  center  of  hemisphere 

2.  polar  point 

c.  Cyclinder  (2  points,  radius) 

1.  center^  of  base£  (2) 

2.  radius  of  cyclinder 

d.  Cone  (2  points,  half-angle) 

1.  apex  point 

2.  center  of  base 

3.  half-angle  (radians) 

e.  Truncated  Cone  (2  points,  half- 
angle,  truncation  point) 

1.  same  as  in  d.  4 

2.  axis  truncation  point 

f.  Hexahedron  (8  points) 

1.  all  apex  points  in  order 
specified  on  diagram 

g.  Ellipsoid  (4  points) 

1.  center  of  ellipsoid 

2.  surface  points  on  3  serai- 
major  axes. 

h.  Truncated  Ellipsoid  (5  or  6  points) 

1.  same  as  in  g. 

2.  axis  truncation  point(s) 

(see  diagram) 

7.  Shield  Type  Code: 

0,1  Hexahedron  2,3  Cyclinder 
4,5  Sphere  6,7  Hemisphere 
8,9  Cone  10,11  Truncated 

12,13  Ellipsoid  Cone 

1 

note:  Even  integers  refer  to  void  shields. 


see  above  for  units  to  be  used. 


Ellipsoid:  Point  Order  and  Location 


Hexahedron:  Point  Order  and  Location 


Figure  23.  Sectoring  Code  Work  Sheet  Instructions 


92 


JMAX  is  further  defined  as  follows: 


JMAX 

Type  of  Ellipsoid  Truncation 

4 

Nontruncated 

+5 

Single  truncation  with  lower  portion 
eliminated  with  respect  to  axis. 

_G 

Single  truncation  with  upper  portion 
eliminated  with  respect  to  +z^  axis 

6 

Double  truncation 

The  shield  types  (ST)  and  number  of  points  required  for  each  shield  type 
are  shown  in  Figure  23. 

The  final  geometric  parameter  is  the  number  of  components  (NCOMP) 
which  comprise  a  composite  shield.  A  composite  shield  represents  a  com¬ 
plex  shield  which  is  defined  by  a  combination  of  as  many  as  ten  solid  and 
void  shields.  Each  void  shield  is  defined  as  voiding  all  the  portions  of  solid 
shields  that  it  intersects.  In  a  composite  shield,  the  all-solid  components 
must  be  followed  by  all  the  void  components  because  of  a  computer  program 
constraint. 

Numerical  coordinates  are  defined  in  a  reference  Cartesian  coordinate 
system  (ABCS)  which  is  fixed  relative  to  the  numerical  model.  The  coordi¬ 
nates  which  define  the  position  and  size  of  the  geometrical  volumes  are  illus¬ 
trated  in  Figure  4.  The  ellipsoid  (x,  y,  z)  axes  must  be  a  right-handed 
orthogonal  set  and  are  not  necessarily  spatially  coincident  with  the  right-hand 
reference  coordinate  system.  Also,  the  z  axis  of  the  ellipsoid  in  the  ABCS 
must  be  the  axis  of  truncation. 

2.  MEVDP  INPUT  DATA  FORMAT 

The  format  for  the  input  data  is  listed  in  Table  VI,  This  table  includes 
the  card  sequence,  subroutines  in  which  the  variables  are  read,  input  varia¬ 
bles,  FORTRAN  nomenclature,  and  the  input  data  format.  Table  VI  also 
includes  the  analysis  and  program  sections  where  the  variables  are  defined 
and  explained  in  greater  detail.  Also,  where  there  is  a  correlation,  ref¬ 
erence  is  made  to  other  input  data  cards.  A  column  designated  as  "Conditions 
for  Data  Inclusion  or  Specific  Parameter  Values"  is  used  to  define  the  con¬ 
ditions  for  which  the  particular  card  is  included  in  the  data  and  the  particular 
values  of  variables  required  for  different  MEVDP  program  computational 
options.  Figure  24  is  a  schematic  showing  the  MEVDP  input  data.  Table  VI 
has  a  reference  column  to  be  used  to  obtain  a  full  explanation  of  the  meaning 
of  a  specific  input  variable. 
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Table  VI 

MEVDP  Input  Data  Format 


I  orm  of  shield  geometrical  configuration 
;  input  data  indit*  or 

a  Card  data  read  and  stored  on 

phvsical  tape  unit  -  10  in  two  files 

b  Data  read  from  reserve  physical 
tape  Unit  10 

c  Card  data  read  and  stored  on 

physical  tape  'Jmt  10  ir.  two  file*  - 
program  execution  terminated 

d  Program  execution  terminated 
|  e.-  Descriptive  title  for  the  geometry  data 
j  ‘  ‘ai  'jr  u«  eript.on  j 


b  Shield  t>pe  (0  -  13,  odd  numbers 
indicate  solid  region,  even  num¬ 
bers  c*'r  "sn  u*  t<  vn»d  g.onl 

c  Material  number  (1  -  50) 

a  Radius  (inches)  of  volume  or 

radian  half-angle  of  cone  (regular 
or  truncated),  whichever  is 
applicable 

e  Material  density  (gm/cm^) 

f  Number  of  sets  of  coordinates  (if 
ST  =  12  or  13.  for  ellipsoid 
1=  4  regula'  ***'»p*«vid  -  h  two 
truncations.  =  +0  upper  portion 
with  respect  to  its  semimajor. 

=  -5  bottom  portion) 

g.  Number  of  elemental  volumes 
defining  the  solid  all  solid  ele-: 
mental  volumes  must  be  listed 
for  the  par*icular  shield  before 
void  regions  are  encountered 

(  ourdma’es  of  volume  in  ABCS' 


a  Hexahed. -on  (ST  =  0  or  1) 

4  cards  Ape*  coordinate 
components 

b.  C>  hnder 

l  card  Centerline  coordi¬ 

nate  comp'  -rnts 

c.  Sphere  (ST  :  4  or  :!] 

1  card  Center  and  one 

point  on  sphere 

d  Hemisphere  iST  =  6  or 

1  card  Center  and  one 

point  on  surface 

e  Cone  <ST  '•  d  or  a) 

1  card  Base  center  and 

apex  point 

f  T runcated  cone  (ST  =  10  or  11) 

2  cards  Base  center  and 

apex,  and  trunca¬ 
tion  point 

g  Ellipsoid  iST  =  12  or  13) 

2  or  Center,  semimajor 

3  cards  axes,  and  trunca-: 

tion  points 

Control  card  which  signals  end  of 
!  spacecraft  geometrical  data 

Astronaut  geometrical  configuration 
NASA-MSC  model  Reference  2) 

a  Repeat  1.2  and  1  3  for  astronaut 
geometry  data 

|  *i  D«  scriptive  cards  which  will  be 
printed  imus*  not  include  anv 


X  (I.  J), 

J  2  1.3. 

I  =  1,  JMAX 


112  MAIN 


>.  (I,  J). 

I  J  =  1.3.  1=  1, 
[number  points 


N10GEN  s  1 


N10GEN  *  0 


I  Figure  24 


Section  IV-3..  Card*  1.2,  1.3,  1.5 


iN  1  U  G  r- 1\  £  1 


File  1  of  Unit  10.  Section  IV- 1. 
Figure  23,  Subsection  111-4.  a. 


J  Figure  4  and  23 


For  order  of  points  see  Figure  1  3 
and  Section  DI-7.  Card  3.  1 


Blank  card  GENTAP 


File  2  of  tape  Unit  10,  Figure  19, 
Subsection  III-4.a. 


in  men  '  ard  sequences  I  2  ard  1.3  a*  -eptated  for  the  remaining  elemental  volumes. 

,s  mar  as  ter  <>f  tnesr  elemental  ge  /metrical  volumes  are  combined  to  form  more  complex  shields,  which  are  designated  as  composite  shields.  The 
>  1  el«  i  .rntal  v  ul  jmt  components  are  used  to  void  the  spatial  region  they  occupy  The  composite  shield  permits  imbedding,  in  which  the  positive 
Id  mav  co-tan  all  or  a  portion  of  a  void  shield  Also,  the  solid  shield  components  may  have  different  material  types  and  densities..  Portions  of 

*  v«*ral  of  the  posit. .e  shield  components  rra>  contain  portions  of  a  single  void  shield.  This  versatility  results  in  a  very  significant  improvement  in  the 

*  onii  tf-  r*  pro«»-ntation  capability  and  reduces  the  magnitude  of  the  shield  synthesis  effort. 


I 

> 


Table  VI  (Continued) 
MEVDP  Input  Data  Format 


Card 

Sequence 

Variable  Deacription 

FORTRAN 

Nomenclature 

Fermat 

Routine  Which 
Reads  Data 

Condition  for  Data 
Inclusion  or  Specific 
Parameter  Value 

References 

(.6 

Control  card  which  aignala  end  of 
aatronaut  geometry 

BLANK 

Blank  card 

GENTAF 

2.  1 

Geometrical  configuration  aelection 
parameter 

OPTN 

112 

MAIN 

Figure  24 

a.  Use  data  on  phyaical  tape  Unit  10 

OPTN  *1 

b.  Uae  new  geometry  data-go  to  Card 
Sequence  1.  1 

OPTN  <1 

2.2 

Detector  position  and  tranaformed 
geometry  parametera 

a.  Detector  coordinatea  in  ABCS 

XDET  (I), 

I  *  1.3 

3E12.8 

TRSHLD 

b.  Spacecraft  geometrical  data  uaage 
indicator 

NS 

110 

Figure  24 

1.,  Spacecraft  data  on  first  file  of 
phyaical  Tape  10  (card  image) 
will  be  uaed. 

NS>0 

2.  Spacecraft  data  will  not  be  uaed. 

NS<0 

c.-  Tranaformed  geometry  output  print 
indicator 

PR  NT 

no 

1.  Transformed  geometry  data 
printed 

PRNT  =  1 

2.,  Transformed  geometry  data 
not  printed 

PRNT  *  0 

2.  3 

Crew  configuration  parameters 

TRSHLD 

a.  Number  of  elemental  volumes  per 
astronaut  (NOMAN  =  34  for  NASA  - 
MSC  model  astronaut) 

NOMAN 

1415 

b..  Number  of  astronauts  in  crew 

NMEN 

c.  Number  of  astronaut  limb  com¬ 
posite  shields  rotated  per  aatro¬ 
naut  (NLIMB(I)  *  1,  for  NASA- 
MSC  model  astronaut) 

NLIMB(I), 

I  *  1.  NMEN 

(1)  NUMB  (I)  =  0 

NMEN  =  0 

(2)  NUMB  (I)  =  1 

NMENfcl 

2.4 

Crew  coordinates  and  orientation 

NMEN>0 

Subsection  IH-4.  a«,  Figure  20 

a.  Astronaut  position  and  orientation 

W(I,J), 

J  =  1.3; 

I  =  1.4 

b.  Astronaut  SMCS  origin  position  in 
ABCS  =  W  (4,  K) 

R4A(K). 

K  *  1,3 

2.5 

Angles  of  rotations  for  crew  limb 
components  (NCOMP  *  8  for  NASA- 
MSC  model  astr*  taut) 

EPSLN(I). 

1*1.,  NCOMP 

NMEN>0 

Figure  19,  Table  V 

Comment 

2. 4  and  2.  5  repeated  for  each 
astronaut 

NMEN  >1 

3.  1 

Hexahedron  point  order 

ORD(I,  J), 

J  *  1,4; 

I  *  1.6 

2413 

ESDOSE 

Figure  13  and  Section  HI-7. 

3  2 

Initiator  for  random  number  generation 
(this  parameter  read  but  not  uaed) 

1R 

112 

3.3 

Special  geometrical,  material,  and  ray 
selection  parameters 

a.  Ray  selection  method  indicator 

RANDOM 

1216 

Section  II- 1 

1..  Systematic  ray  selection 

RANDOM  *  0 

2.  Random  ray  selection 

RANDOM  >0 

b.  Number  of  rays  selected  (S1000) 

N 

c.  Number  of  material  types  (350) 

NMAT 

d.  Number  of  shield  elemental  volumes 
requiring  new  densities  relatlva  to 
data  or  physical  tape  Unit  10. 

NRHO 

e.  Direction  cosines  print  indicator 

NPRT1 

1 Direction  cosine  data  printed 

NPRTl  •  1 

2,  Direction  cosine  data  not 
printed 

_ 

NPRT1*  0 
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Table  VI  (Continued) 
MEVDP  Input  ^ata  Format 


Ca  rd 
Sequence 

Variable  Description 

FORTRAN 

Nomenclature 

Format 

Routine  Which 
Reads  Data 

Condition  for  Data 
Inclusion  or  Specific 
Parameter  Value 

References 

f.  Detail  ray  tracing  print  indicator 

NPRT2 

1.  Tracking  data  printed 

NPRT2  =  1 

Section  i'-l 

2.  Tracking  data  not  printed 

nprt?  =  o 

g  Number  of  macroscopic  solid-angle 
increments  for  systematic  '■ay 
aelectior  (5l0) 

NDIV 

h.  Sorted  material  and  path-length  out¬ 
put  v*orage  medium  indicator 

IPUNCH 

Sec*  ion  IiI-19  and  IV- 3. 

1.  Sorted  data  printed  only 

IPUNCH  =  0 

2.  Sorted  data  stored  on  tape  Unit  7 
with  format  for  punched  card 
output  and  printed 

PUNCH  =  7 

3.  Sorter  data  on  tape  Unit  8  m 
binary  form 

IPUNCH  =  8 

4.  Sorted  data  stored  on  Units  7 
(punched  card  format)  and  8 
(binary  form)  and  printed  1 

0<  IPUNCH <7 

3  4 

Shield  serial  numbers  for  shields  with 
specified  density  changes 

NSS(I). 

I  =  1.  NRHO 

1216 

NRHO  >0 

Card  3.  3d 

3  5 

New  shield  densities  corresponding  to 
shield  serial  numbers 

DENS(I), 

1=1.  NRHO 

6E12.8 

NRHO>0 

Card  3.  4 

3.6 

Title  card  for  problem  description 

CASE 

1 8A4 

ESDOSE 

3.  7 

Curve  fit  data  for  standard-material 
areal-density  curve  fit  parameter 

a  Tape  9  usage  parameter 

NEWF9 

112 

Section  IV- 3. 

l.'  Store  curve  fit  data  on  new 
reserve  tape  on  tape  Unit  9 
for  one  or  several  dosimeter 
positions 

NEWF9  =  -I 

2.  Curve  fit  data  will  not  be 
saved  on  tape 

NEWF9  =  0 

3.  Curve  fit  data  added  to  a 
previously  generated  tape 
which  is  mounted  on  tape 

Unit  9.  If  the  data  are  to  be 
added  after  File  I  for  several 
dosimeter  positions,  then 

NEWF9  *  I,  1  ♦  1.,  I  +  2,  etc. . 
for  the  1,2,3.  etc..,  dosimeter 
locations. 

NEWF9  >0 

b.  Cut  off  equivalent  areal  density  for 
ray  tracing,  including  sorted  data 

’MX 

£12.8 

c.  Parameter  used  to  compensate 
for  tracking  roundoff  error  (sug¬ 
gested  value  =  0.001) 

EPSLN 

£12.  8 

3.  8 

Range-energy  curve  fit  parameters 

6E12.8 

Subsection  11*7.  a. 

a.  Standard-material  parameters 
(6STD.  tSTDl 

DSTD,  ESTD 

b.  Actual  material  parameters 
(6t,  V  i  =  1,,  NMAT) 

DELTA(I), 

ETA  (I), 

1=1.,  NMAT 

3.  9 

Standard-material  name  for  print 
output 

TTLSTD(I), 

I  =  1.3) 

3A10 

3.  10 

Actual  material  names  for  print 
output  (two  per  card  with  30  col¬ 
umns  each) 

TTLMAT  (I. 

J).  1=1.3. 

J  =  1 ,  NMAT 

6A10 

3.  1 1 

Systematic  ray  selection  parameters 

ESDOSE 

RANDOM  =  0 

Subeection  11- 1 .  b. 

a  Number  of  microscopic  solid- 
angle  discrete  divisions 

N 

112 

b  Azimuthal  angles  (deg.  ees) 

PHII,  PHIF 

4E12.  8 

c  Colatitudv  angles  0j,  0p  (degrees) 

THETI. 

THETF 

Comment 

Repeat  3  1  1  for  NDIV  macroscopic 
solid  angles  with  Z.  Nc  s  1000  where 
•\c  «  number  of  solid  angles  computed 
iby  the  program  corresponding  to  N 

3.  12 

Areal-density  curve  fit  sections 
parameters 

CEOMDS 

Section  1Q«18. 

a.  Number  of  sections  (55) 

NSECT 
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Table  VI  (Concluded) 
MEVDP  Input  Data  Format 


Card 

Sequence 

Variable  Description* 

FORTRAN 

Nomenclature 

Format 

Routine  Which 
Reads  Data 

Condition  for  Data 
Inclusion  or  Specific 
Parameter  Value 

Refe  rences 

b.  Degree  of  polynomial  used  for  each 
flection  (5 19  per  section) 

NSDEG(l), 

1=1.  NSECT 

6112 

3.  13 

Areal-density  curve  fit  point  density 
parameters 

4E12.8 

GEOMDS 

Section  III-  18 

a.  Maximum  areal  density 

XMAX  TMX 
(3.7b) 

b.  t  (suggested  value  *  8.0) 

TAU 

c.  iR  (suggested  value  -  0.  105) 

DELR 

d.  Maximum  areal  density 

TMAX  = 

XMAX 

3.  14 

End  of  data  case  for  a  dosimeter 
position  (IDUMMY  must  be  less 
than  zero) 

IDUMMY 

112 

Comment 

GO  to  Card  2.  1 
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TITLE  CARD 


i 


L 
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Figure  24.  MEVDP  Input  Data  Schematic 


3.  TAPE  UNIT  UTILIZATION 


Table  VII  is  a  summary  of  the  tapes  required  in  the  MEVDP;  individual 
and  more  detailed  explanation  of  the  data  on  each  tape  follows. 

Table  VII 

Tape  Unit  Use  and  Requirements 


Tape  Unit 

Use 

5 

Input 

6 

Output 

7 

Storage  •*  punch  card 
output  of  sorted  data 

8 

Storage  -  binary  form 
of  data  on  Tape  7 

9# 

Multiple  file  storage  tape  - 
optional  physical  tape 

10 

Input  data  -  optional  physical 
tape  -  (See  cards  1.  2  and  1.  3, 

Table  VI) 

12 

Scratch  disk  file 

13 

Scratch  disk  file 

*Note: 

Tape  9,  a  physical  tape,  is  only  required 
for  permanent  storage  of  the  standard- 
material  areal-density  distribution  data. 

The  computer  system  tape  units  5  and  6  are  used  as  standard  input  and 
output  devices.  For  the  other  required  tape  units,  the  following  definitions 
will  be  used: 

SYSDA  =  system  direct  access  device  (disk  scratch  unit) 

Temporary  storage  unit  =  disk  unit 

Physical  tape  unit  =  actual  mounted  reserve  or  scratch  physical 

tape  unit,  a  physical  unit  tape 

Unit  10  contains  two  card-image  files:  (1)  spacecraft  geometry  and 
material  model,  and  (2)  the  astronaut  geometry  and  material  model.  The 
numerical  geometrical  model  can  be  input  as  cards  and  written  on  disk 
unit  10,  or  unit  10  can  be  a  mounted  reserve  tape  containing  the  required 
two  numerical  model  files.  For  additional  discussion,  refer  to  cards  1.  1 
through  2.  2  of  Table  VI,  and  Subroutine  GFNTAP  (Section  III-2.  ). 


99 


Unit  12,  a  scratch  disk  file,  is  used  for  temporary  storage  of  numerical 
model  data  after  they  have  been  transformed  to  the  RDCS  for  an  input 
dosimeter  position. 

Unit  13,  a  scratch  disk  file,  is  used  for  temporary  storage  of  the 
astronaut  numerical  model  data  (NASA  -  MSC  Astronaut  Model)  after  they  have 
been  transformed  to  the  ABCS.  Subroutine  TRSHLD  uses  unit  13  to  generate 
the  transformed  data  for  the  crew  in  the  RDCS,  and  adds  the  data  to  the 
transformed  spacecraft  data  on  unit  12.  Other  references  to  unit  13  are  in 
Subroutine  TRSHLD  (Section  III-4.  ). 


Unit  9,  a  physical  unit  file,  is  used  for  storage  of  the  standard-material 
real  density  distribution  curve  fit  data.  If  the  curve  fit  data  are  not 
requested,  unit  9  can  be  specified  as  a  scratch  disk  file.  In  general,  unit  9 
will  be  a  multifile  reserve  tape  with  following  data  records  for  each 
dosimeter  position: 

Record  1. 

a. 

NRAYS 

=  number  of  rays  for  tracking 

b. 

NPTS 

=  number  of  points  used  for  the 
curve  fits 

c. 

NSECT 

=  number  of  curve  fit  sections 

d. 

[NDEG  (I),  I  = 

1,  NSECT] 

=  degree  of  least-square  polynominal 
used  for  each  section 

Record  2. 

a. 

TMIN 

=  minimum  equivalent  thickness 
(grams/cm")  computed  for  a 
particular  ray 

b. 

TMAX 

=  specified  maximum  equivalent 
thickness  for  ray  tracing 

c. 

[XFIT(I),  I  =  1,  NSECT] 

=  maximum  equivalent  thickness  for 
each  sectional  curve  fit 

d. 

[THK(J),  J  =  1 

,  NPTS] 

=  computed  weighted  equivalent 
thickness  values  for  the  areal- 

density  distribution  function 


100 


Record  3. 


[PATH(I),  1=1,  NRAYS]  =  total  equivalent  thickness  for  each  ray 

Next  NSECT  Records 

The  generalized  nth  record  contains  the  following  polynominal  curve 
fit  coefficients: 

[CSECT(J,N),  J  =  1,  NDEG(N)  +  1]  =  polynominal  curve  fit  coefficients 

For  further  discussion  refer  to  card  3.  7  of  Table  VI  and  subroutine 
GEOMDS  (Section  3.  18). 

The  areal-density  function  fractional  solid  angle  (Y^)  corresponding  to 
the  ith  areal-density  point  value  T^  is  computed  from  the  curve  fit 
coefficients  with 


NDEG(N)  +  1 

Yt  =  ^  CSECT  (R,N1  Tj*"1 

R  =  1 


where 

XFIT(N-l)  <  Tt  <  SFIT(N) 

Unit  7,  system  punch  card  output  tape,  and  unit  8,  a  scratch  disk  file 
or  a  physical  reserve  tape,  are  used  for  storage  of  the  sorted  computed  path- 
length  and  material  data.  These  output  units  are  used  according  to  the  value 
of  the  variable  IPUNCH  as  given  in  Table  III. 


Table  VIII 

Use  of  Output  Tape  Units  Versus  IPUNCH 


IPUNCH 

So 

rted  Data  Output  Form 

Print  Data 

On  Unit  6 

Punch  Data 

On  Unit  7 

Binary  Data 

On  Unit  8 

0 

Yes 

No 

No 

7 

Yes 

Yes 

No 

8 

Yes 

No 

Yes 

0<  IPUNCH<  7 

Yes 

Yes 

Yes 
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For  tape  unit  7,  the  IRth  record  of  th**  first  NRAYS  records  contain  the 
macroscopic  solid-angle  parameters. 

IK.  [DIRCOS(J.IR),  J  =  1,  3],  NRSA,  FSA 


with  the  formal 


15,  3E15.  8,  15,  El 5.  8 


whe  re 


IR  -  ray  number 

DIRCOS  =  ray  direction  cosine  components 

NRSA  =  number  of  rays  for  the  macroscopic  solid  angle  containing 
the  IRth  ray 

FSA  =  solid  angle  subtended  by  the  IRth  ray 

These  records  are  followed  by  the  sorted  material  and  path-length  data. 
Each  record  contains  the  variables 

IR,  I,  MAT(I),  X(I),  DX(I) 


with  the  format 


315,  2E15.  8 


where 


IR  =  ray  number 
I  =  path-length  number 
MAT (I)  -  material  number 

X(I)  =  sorted  ordered  distance  for  the  Ith  path  of  the  IRth  ray 

DX(I)  =  traversal  areal  density  for  the  Ith  path  length  of  the  IRth  ray 

The  macroscopic  solid-angle  parameters  and  sorted  path-length  and 
material  data,  which  can  be  written  on  unit  7,can  also  be  written  on  unit  8  in 
binary  form,  i.  e.  ,  Without  format.  Arrangement  of  the  data,  by  record,  is 
the  same  for  units  7  and  8.  Additional  discourse  on  the  use  of  units  7  and  8 
is  in  card  3.  3h  of  Table  VI,  and  subroutine  ORDER  (Section  III- 19). 
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SECTION  V 


SAMPLE  PROBLEM  SOLUTION 


A  sample  problem,  supplied  by  AFWL,  is  defined  in  Figures  25  and 
26.  The  elemental  volume  components  of  the  sample  problem  are  depicted 
in  Figure  25  and  the  composite  configuration  is  illustrated  in  Figure  26  with 
associated  dimensional  parameters.  Tables  XI  and  XII  tabulate  (1)  the 
material  properties,  (2)  geometrical  dimensional  data,  and  (3)  the  derived 
coordinates  and  parameters  required  as  input  data  for  the  MEVDP.  The 
dimensional  parameters  in  Table  X  are  illustrated  in  Figure  26. 

Input  data  for  the  sample  problem  are  listed  in  Table  XII  for  three 
dosimeter  locations.  The  cards  are  referenced  to  the  input  data  format  in 
Table  VI  and  the  input  data  schematic  in  Figure  24  through  the  identification 
field,  i.  e.  ,  columns  73  to  80.  The  numbers  in  columns  73  through  76  refer 
to  the  card  numbers,  for  example  31  and  314  are  interpreted  as  cards  3.  1 
and  3.  14,  respectively.  Lines  1  through  53  also  indicate  the  corresponding 
shield  serial  numbers  for  the  sample  problem  geometrical  configuration. 
Table  IX  shows  the  correlation  between  the  shield  serial  numbers  and  the 
shield  elemental  volume  number  in  Figure  25. 

In  order  to  assure  the  accuracy  of  the  transformation  of  the  AFWL 
sample  problem  data  into  proper  format  for  the  MEVDP,  cathode  ray  tube 
(CRT)  pictures  were  generated  for  several  cross  sections.  These  cross  sec¬ 
tions  are  depicted  in  Figures  27  through  30.  The  titles  denote  the  planes  of 
the  cross  sections  and  their  center-point  coordinates  in  inches  relative  to  the 
coordinate  system  in  Figure  25.  The  sample  problem  data  have  be^r  execu¬ 
ted  with  similar  results  on  the  North  American  Rockwell  Space  Division 
(NR  SD)  IBM  System  360  and  the  AFWL  CDC  System  6600.  The  insignificant 
numerical  differences  were  due  to  inherent  differences  in  the  number  of  sig¬ 
nificant  figures  for  the  two  computer  systems. 

This  section  is  concluded  with  a  brief  explanation  of  the  output  data 
format.  Appendix  lisa  listing  of  excerpts  from  the  NR  SD  IBM  System  360 
computer  solution  for  the  AFWL  sample  problem  with  the  dosimeter  at 
(x  =  0,  y  =  0,  z  =  0  inches). 

There  were  two  changes  in  the  program  listing  (Appendix  II)  in  sub¬ 
routine  ORDER.  The  dimension  of  BPL  (card  F2)  and  the  variable  NMAX 
(card  F12)  were  reduced  to  3000,  This  change  was  made  to  demonstrate  the 
capability  of  the  program  to  perform  an  internal  sort  of  the  path-length  data 
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POINT 


Elemental  Volume  Components  of  the  AFWL  Sample  Problem 


26.  AFWL  Sample  Problem  Shield  Configuration 


Table  IX 


Table  X 

AFWL  Sample  Probleni  Dimensions  and  Dosimeter  Locations 
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Table  XI 

Shield  Coordinates  and  Parameters 


Table  XI  (Concluded) 

Shield  Coordinates  and  Parameters 
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0001  1  AFHL  SAMPLE  PROBLEM  FOR  THE  MEVDP 
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COLUMN  NUMBERS 


Table  XII  (Continued) 

AFWL  Sample  Problem  Input.  Data 
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Table  XU  (Continued) 

AFWL  Sample  Problem  Input  Data 
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Table  XII  (Concluded) 
AFWL  Sample  Problem  Input  Data 
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Figure  27.  AFWL  Sample  Problem  (x,  y)  Plane,  Center  (0,  0,  0) 
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Figure  28.  AFWL  Sample  Problem  (y,  z)  Plane,  Center  (0,  0,  0) 


Figure  29.  AFWL  Sample  Problem  (y,  z)  Plane,  Center  (-20. 1, 0, 0) 


when  the  amount  of  data  exceeds  the  available  core  space.  The  variable 
NMAX  is  three  times  the  number  of  path  lengths  for  which  the  associated 
data  can  be  sorted  in  the  available  core  space. 

The  output  data  commence  with  a  descriptive  title  for  the  input  geo¬ 
metrical  configuration  data  which  follow.  The  geometrical  data  include  the 
record  number  and  associated  card  image  of  the  input  geometry  data  which 
are  stored  on  physical  tape  unit  10.  These  data  correspond  to  the  data  for¬ 
mat  of  cards  1.  2  and  1.  3  in  Table  VI. 

The  geometrical  configuration  data  are  followed  by  the  card  images 
for  the  astronaut  geometry  data.  Since  NMEN  =  0  (see  card  1.  5b  in 
Table  VI),  there  are  no  astronauts,  and  the  astronaut  geometrical  records 
only  include  a  descriptive  ootput  message. 

The  next  output  data  are  the  dosimeter  coordinates  in  the  ABCS. 
Because  the  variable  PRNT  is  unity,  ‘■he  transformed  geometry  data  are 
printed  next  (see  card  2.2  of  Tables  v'l  and  IX).  These  data  include  (1)  the 
shield  serial  number,  (2)  shield  geometrical  form  or  type,  (3)  shield  octant 
location  in  the  RDCS,  (4)  material  type,  and  (5)  the  material  density  in 
grams  per  cubic  centimeter.  The  transformed  data  also  include  the  coor¬ 
dinates  (inches)  of  the  defining  points  of  each  shield  in  the  RDCS  and  the 
rotation  matrix  which  transformed  the  point  coordinates  from  the  DSCS  to 
the  RDCS.  Other  miscellaneous  data  include  the  sphere  and  cone  radii, 
ellipsoid  major  axes,  the  cone  half-angle  in  radians  and  the  cone  base  (RB) 
and  truncation  radius  (RT).  These  coordinates  and  length  dimensions  have 
inch  units. 

The  transformed  data  are  followed  by  a  listing  of  the  number  of  shields 
in  each  octant  in  the  RDCS,  and  the  CPU  time  expended  for  data  transforma¬ 
tion  in  minutes. 

A  descriptive  problem  title  is  printed  next,  followed  by  the  macro¬ 
scopic  solid  angle  «t>  (PHI)  and  6  (THETA)  boundaries  in  degree  units  in  the 
DSCS.  The  next  section  is  a  listing  of  the  ray  number,  colatitude  and 
azimuthal  angles  of  each  ray  direction  in  the  DSCS,  and  the  direction  cosines. 
These  data  are  not  printed  if  NPRTI  =  0  (see  card  3.  3  of  Table  VI). 

The  ray  directional  data  are  followed  by  the  number  of  rays  in  the 
particular  macroscopic  solid-angle  increment  and  the  subtended  solid  angle 
per  ray  in  steradians. 

The  next  output,  comprising  three  columns,  denotes  the  material 
range-energy  curve  fit  parameters  (see  Equation  49)  6  (DELTA)  and  T  (ETA), 
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and  the  material  names.  The  abbreviation  "STD"  refers  to  the  standard 
material,  and  the  subscripted  variables  refer  to  the  actual  materials  which 
are  used  in  the  geometrical  configuration  data  cards  for  the  material  types. 
The  subscript  I  refers  to  the  Ith  material  type. 

The  material  data  are  followed  by  detail  path-length  data  for  each  com¬ 
posite  shield  if  NPRT2  =  1  (see  card  3.  3  of  Table  VI).  These  data  include 
the  number  of  composite  shield  components,  and  the  serial  number,  geomet¬ 
rical  type,  and  density  (grams  per  cubic  centimeter)  for  each  component; 
the  ray  numbers,  direction  cosines  in  the  DSCS,  path-length  areal  density 
(XS  in  grams  per  square  centimeter),  and  the  accumulative  areal  density  for 
the  ray  directions.  These  output  data  occur  only  for  rays  which  intercept  the 
composite  shield. 

The  "Track  Time"  is  the  CPU  time  in  minutes  used  to  ray  trace  the 
elemental  volumes. 

Next,  the  total  areal  density  is  printed  for  the  standard  material  for 
each  ray  direction.  Each  row  starts  with  the  areal  density  for  the  ray  number 
in  the  left  column  and  corresponds  to  increasing  ray  numbers  moving  from 
left  to  right. 

The  ray  areal  density  is  followed  by  the  areal-density  curve  fit  para¬ 
meters,  in  Equation  60,  i.  e. 


XMAX  =  T  ,  AV 
MAX 

TAU  =  t 
DELR  =  6R 

These  are  in  turn  followed  by  the  weighted  values  of  thickness  (areal  density 
in  grams  per  square  centimeter)  computed  with  Equation  61 . 

The  next  set  of  output  data  is  the  computed  fraction  of  solid  angle 
("Distribution"  column)  with  areal  densities  less  than  the  corresponding  areal 
densities  ("Thickness"  column).  The  fractions  of  solid-angle  values  are  com¬ 
puted  from  the  standard-material  areal  densities  for  each  ray  and  the  weighted 
values  of  areal  density.  The  column  entitled,  "Curve  Fit  Points",  is  the 
curve  fit  of  fraction  of  solid  angle  i.  e. ,  the  "Distribution"  column.  The 
least-square  curve  fit  is  computed  with  Equation  59.  The  subsequent  three 
columns  of  data  entitled,  "Coefficients"  are  the  values  of  the  Cj  constants  in 
Equation  59  where  the  values  of  j  increase  row-wise  from  left  to  right.  This 
type  of  data  is  repeated  for  each  of  the  specified  number  of  curve  fit  sections 
(see  card  3. 12  in  Table  VI),  i.  e. ,  three  for  this  case. 
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The  CPU  times  used  in  subroutines  GEOMDS  and  ESDOSE  are  printed 

next. 

The  output  data  are  concluded  with  listings  of  the  unsorted  and  sorted 
data  for  each  path  length  of  material  encountered  for  each  ray.  The  next 
line  of  output  data  denotes  the  total  number  (1604)  of  path  lengths  encountered 
for  all  the  rays  and  the  number  (998)  of  path  lengths  which  can  be  sorted  with 
the  available  core  space,  which  normally  corresponds  to  7,  000  path  lengths 
and  was  changed  to  1,  000  for  this  problem. 

The  next  set  of  data  is  a  listing  of  the  unsorted  data,  denoting  the  (1) 
ray  number,  (2)  path  number  for  the  ray,  (3)  distance  to  the  path  length, 

(4)  path  length  or  distance  through,  and  (5)  the  path-length  material  number. 
The  unsorted  data  are  followed  by  a  listing  of  the  areal  densities  sorted 
according  to  ray  number  and  distance  to  the  path  lengths  traversed.  The 
sorted  data  include  the  (1)  ray  number,  (2)  path-length  number  for  the  ray, 

(3)  distance  to  the  path  length,  (4)  path-length  areal  density,  (5)  path-length 
material  number,  and  (6)  the  material  name.  The  line  in  the  middle  of  the 
sorted  data,  denoting  606  of  1604  path  lengths,  indicates  that  the  unsorted 
data  tape  had  to  be  read  a  second  time  because  of  insufficient  core  space 
allocation  for  sorting  all  the  path  lengths.  The  sorted  data  follow  for  the 
remaining  606  path  lengths. 
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APPENDIX  I 

MEVDP  AFWL  SAMPLE  PROBLEM 
SOLUTION  OUTPUT  DATA 
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AFWL  SAMPLE  PROBLEM  FOR  THE  MEVOP 
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TRACKING  OF  PROTONS  THROUGH  A  HEMISPHERE 
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Figure  18.  Functional  Diagram  of  the  MEVDP 


